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AERONAUTIC SYMBOLS 
1. FUNDAMENTAL AND DERIVED UNITS 



Symbol 

Metric 

English 

Unit 

Abbrevia- 

tion 

' 

Unit 

Abbreviation 

Length 

Time 

l 

t 

F 

meter.. A — 

second. ... 

m 

8 

kg 

foot (or mile) 
second (or hour) 

ft (or mi) 
sec (or hr) 
lb 

Foroe 

weight of 1 kilogram 

weight of 1 pound. 


P 

V 

horsepower (metric) 


horsepower. 

hp , 

mph 

fps 

Speed 

/kilometers per hour 

\meters per second 

kph 

mps 

milea per hour 

feet per second.. . . 


2. GENERAL SYMBOLS 


Weight =m <7 

Standard acceleration of gravity=9.80665 m/s 2 
or 32.1740 ft/sec 2 
W 

Mass= — 

9 

Moment of inertia =mfc 2 . (Indicate axis of 
radius of gyration k by proper subscript.) 
Coefficient of viscosity 


v Kinematic viscosity 

p Density (mass per unit volume) 

Standard density of dry air, 0.12497 kg-m~ 4 -s 2 at 15 c 
and 760 mm; or 0.002378 lb-ft~ 4 sec 2 
Specific weight of “standard” air, 1.2255 kg/m 8 
0.07651 lb/cuft 


C 


or 


3. AERODYNAMIC SYMBOLS 


Area 

Area of wing 
Gap 
Span 
Chord 

Aspect ratio, ^ 

True air speed 
Dynamic pressure, ^ pV 1 

Lift, absolute coefficient 

Drag, absolute coefficient Co=|g 

Profile drag, absolute coefficient C Dq 


Do 

qS 

D t 


Induced drag, absolute coefficient C Di =-^g 


D v 


Parasite drag, absolute coefficient C Dp —-^ 


i„ Angle of setting of wings (relative to thrust line) 
it Angle of stabilizer setting (relative to thrust 
line) 

Q Resultant moment 

Q Resultant angular velocity 

VI 

R Reynolds number, p — where l is a linear dimen- 
sion (e.g., for an airfoil of 1.0 ft chord, 100 
mph, standard pressure at 15° C, the corre- 
sponding Reynolds number is 935,400; or for 
an airfoil of 1.0 m chord, 100 mps, the corre- 
sponding Reynolds number is 6,865,000) 
a Angle of attack 

e Angle of downwash 

«o Angle of attack, infinite aspect ratio 

a< Angle of attack, induced 

a a Angle of attack, absolute (measured from zero- 

lift position) 

y Flight-path angle 


G 


Cross-wind force, absolute coefficient 
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A RECURRENCE MATRIX SOLUTION FOR THE DYNAMIC RESPONSE OF AIRCRAFT IN GUSTS 1 

By John C. Houbolt 


SUMMARY 

A systematic 'procedure is developed jor the calculation of the 
structural response of aircraft flying through a gust by use of 
difference equations and matrix notation. The use of differ- 
ence equations in the solution of dynamic problems is first 
illustrated by means of a simple-damped-oscillator example. 
A detailed analysis is then given which leads to a recurrence 
matrix equation for the determination of the response of an 
airplane in a gust. The method takes into account wing 
bending and twisting deformations, fuselage deflection, vertical 
and pitching motion of the airplane, and some tail forces. The 
method is based on aerodynamic strip theory, but compressi- 
bility and three-dimensional aerodynamic effects can be taken 
into account approximately by means of over-all corrections. 
Either a sharp-edge gust or a gust of arbitrary shape in the 
spanwise or flight directions may be treated. In order to aid in 
the application of the method to any specific case, a suggested 
computational procedure is included. 

The possibilities of applying the method to a variety of tran- 
sient aircraft problems, such as landing, are brought out. A 
brief review of matrix algebra, covering the extent to which it is 
used in the analysis, is also included. 

INTRODUCTION 

In the problem of an airplane flying through gusts, accurate 
predictions of stresses are not always obtained if the inter- 
action between aerodynamic loads and structural deforma- 
tions is not considered. The present report gives a method 
for determining the dynamic response of aircraft in gusts in 
which this interaction is considered. An approach is em- 
ployed which is a departure from the usual modal type of 
solution. The time derivatives in the integrodifferential 
equations of motion of the airplane are replaced by appro- 
priate difference expressions and use is made of matrix nota- 
tion to express conveniently the conditions of equilibrium at 
a number of points along the wing span. The result is a 
systematic procedure which is complete and general in form. 
The airplane is assumed to be free to translate and pitch. 
Wing bending, wing twist, and fuselage flexibility are all in- 
cluded. Tail forces due to vertical motion, angle of attack, 
and gust penetration are also included in the analysis. 

With the method, a gust with any gradient in the direction 
of flight or along the span may be treated. The method is 


based on aerodynamic strip theory, but over-all compressi- 
bility and aspect-ratio corrections may be included with- 
out difficulty, if desired. One such over-all correction is 
indicated. 

In the first part of the report the method of using difference 
equations in the solution of dynamic problems is illustrated 
by an example in which the transient response of a simple 
oscillator is determined. The analysis for the determination 
of the response of an airplane in a gust is then given. In the 
following section a computational procedure is suggested. 
This section is not intended to describe or add to the under- 
standing of the analysis, but by following the directions indi- 
cated, the response of any airplane may be found without 
following through the complete details of the analysis. 

Supplementary definitions and derivations are presented in 
appendixes. Appendix A summarizes the various matrix 
coefficients and matrices that are used in the analysis, 
appendix B gives a derivation of the difference equations, 
appendix C gives a derivation of the flexibility matrices, 
appendix D gives a derivation of a recurrence equation for 
evaluating the form of Duhamel’s integral which involves an 
exponential kernel, and appendix E presents a review of the 
fundamentals of matrix algebra. It is suggested that those 
not familiar with matrix algebra read appendix E before 
reading the analysis. 

SYMBOLS 

a distance between leading edge of wing and elastic 

axis 

a i coefficient used in unsteady lift function for sudden 

change in angle of attack 
A aspect ratio of wing 

A t aspect ratio of horizontal tail 

b semispan of wing 

c chord of wing 

c 0 chord at wing midspan 

c 0t midspan chord of tail 

c t mean aerodynamic chord of tail 

e distance between mass center of wing cross section 

and elastic axis of wing, positive when elastic 
axis lies forward of mass center 
E Young’s modulus of elasticity 


1 Supersedes NACA TN 2060, “A Recurrence Matrix Solution for the Dynamic Response of Aircraft in Gusts’' by John C. Houbolt, 1950. 
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suddenly applied force 
shear modulus of elasticity 

integers 0, 1, 2, 3, 4, and 5 used to designate sta- 
tions (for most part used as parenthetical num- 
bers, that is, w(3) is deflection at station 3) 
bending moment of inertia 
torsional stiffness constant 

radius of gyration of wing mass about elastic axis 
or elastic spring constant 

length of section associated with a spanwise station 
aerodynamic lift over interval l on wing 
shear force transmitted to wing by fuselage 
aerodynamic lift over interval l on whig due to gust 
one-half aerodynamic lift on tail due to gust 
one-half total aerodynamic lift on tail 
part of aerodynamic lift over interval l on wing (see 
equation (16)) 

part of aerodynamic lift over interval l on wing (see 
equation (17)) 

mass of beam included in interval l or concentrated 
mass in spring oscillator 

mass m including apparent mass effect (m + — 

assumed over-all compressibility and aspect-ratio 

correction for wing ( ^ 

\2 + ylyl-M 2 / 

assumed over-all compressibility and aspect-ratio 
correction for horizontal tail ( 

\2 + A t y'l-M*J 

mass moment me including apparent mass effect 

/ , tt pic 3 ( 1 a\\ 

(me+— 

mass of fuselage per unit length 

mass polar moment of inertia mk 2 including appar- 


(HT 


ent mass effects ^mk 2 + wp ^ c (i ^~ ) 

Mach number or aerodynamic moment over inter- 
val l about elastic axis of wing 
moment transmitted to wing by fuselage 
integers 0, 1, 2, 3, and so forth to designate num- 
ber of time intervals passed 
normal load acting at a station 
fuselage inertia loading per unit length 
torsional load acting at a station 

distance traveled by wing in half-chords ( “ — t, 

\ Co 

where midspan chord c 0 is used as reference 
cliord^ 

distance interval in half-chords 
horizontal-tail area 

time, zero at beginning of gust penetration 
forward velocity of flight 
vertical velocity of gust 


deflection of elastic axis of wing, positive upward, 
or deflection of mass oscillator 
deflection of fuselage, positive upward 
fuselage modal function, zero at wing elastic axis 
and unity at tail one-quarter-chord location 
distance along fuselage measured from wing elastic 
axis, positive in rearward direction 
distance from foremost part of nose to elastic axis 
distance from elastic axis to one-quarter-chord loca- 
tion on tail 

distance along wing measured from center of air- 
plane 

ratio of dynamic deflection to static deflection 
angle of attack of horizontal tail, positive in the 
stalling direction 

forward-speed and aspect-ratio factor for wing 
(: m A irpU ) or coefficient of damping for spring 
oscillator 

forward-speed and aspect-ratio factor for tail 

Q m At irpS, U^J 

exponential coefficient in 4> function associated with 

/ 2(7 \ 

time t (r= — 

coefficient of fuselage modal function 
time interval 

exponential coefficient in $ function associated 
with variable s 

dimensionless interval between i — 1 and i stations 
(' \ib is actual length) 
mass density of air 

angle of twist of wing, positive in stalling direction 
function which denotes growth of lift on rigid air- 
foil entering sharp-edge gust (used without sub- 
script to indicate function for wing and with 
subscript t used to indicate function for tail) 
natural frequency associated with IV U radians per 
second 

unit-step function 

function which denotes growth of lift on airfoil 
following sudden change in angle of attack (used 
without subscript to indicate function for wing and 
with subscript t used to indicate function for 
tail) 

square matrix 
rectangular matrix 
column matrix 
row matrix 


Subscripts: 

t 

0, 1, 2, 3, . . . n 
0, 1, 2, 3, 4, 5 or i 


number of time intervals passed 
station (however, station is usually 
given as parenthetical number, such 
as w(3) for deflection at station 3); 
i is also used as general subscript in 
appendix A 
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All the terms, coefficients, and matrices not defined in this 
section are defined in appendix A. 

Dots are used to indicate derivatives with respect to time; 

. , bw . bw 

for example, =w or ^~ =w ' 

ANALYSIS 

TRANSIENT RESPONSE OF A SIMPLE DAMPED OSCILLATOR 

In order to illustrate the use of difference equations and 
to test the accuracy of the procedure that is to be used in 
the solution of the more complicated gust problems, the 
solution of a simple problem having a known analytical 
solution is first presented. The problem is to compute the 
response of the damped oscillator shown in figure 1 to a 
suddenly applied force. The differential equation of motion 
of this system due to the suddenly applied force is 


m w + /3w + kw = FI (t) 


( 1 ) 


By use of difference equations this differential equation may 
be transformed into an equation which involves deflection 
ordinates at several successive values of time. Probably 
the most commonly used difference equations are the follow- 
ing (see appendix B for derivation): 


Wn 


W„ + l— W„ 


2e 

w„+, — 2w„+w„_, 


( 2 ) 


(3) 


which give the derivatives at the intermediate of three 
successive ordinates. Although these equations are quite 
adequate for the oscillator problem of the present report, 
they cannot be used in the gust analysis which follows. 
Rather, for reasons which are brought out in a subsequent 
part of the analysis, equations that give the derivatives at 
the end ordinate of several successive ordinates must be 
used. If only three successive ordinates are used, the deriv- 
atives so found are not accurate enough to be useful. If 
a fourth ordinate is added, however, derivatives may be 
taken at the end ordinate with accuracies which are com- 
parable to those given by equations (2) and (3). Such 
derivatives are derived also in appendix B and are given 
by the equations: 




W n 


l\w n — 1.8w n _i + 9w„_2— 2 w„- 3 
6e 

2w„ — 5w n _,-f 4w n _ 2 — 


(4) 


(5) 


Although either equations (2) and (3) or equations (4) and 
(5) may be used in the solution of this oscillator problem, 
only equations (4) and (5) will be used, since only these 
equations can be used in the gust-problem solution presented 
in this report. 

If the derivatives in equation (1) are replaced by the 



F Kin re 1. — Damped, oscillator and suddenly applied force. 

difference equations (4) and (5), the following equation is 
obtained: 

( 2+ T m + ^ e 2 ) W “ K 5 + ^) W ' , - 1_ ( 4 + ^) W,! - 2 + 

w„ 


( 


1+MC 

T 6m, 


3 + 


IF 

m 


( 6 ) 


This equation may be said to be the difference equation of 
motion. If the more general case of a variable applied force 
were being considered, the factor F in this equation would 
be replaced by F v , the value of the force present at the 
time t—nt. 

k 

If a specific case is now considered, in which —=400, 

-^- = 2, 6 = 0.01, F= 1, and the notation z — ^wrr (ratio of 
2 m r/lc 

dynamic deflection to static deflection) is used, equation (6) 

becomes 


z,r- 


= 0.0 18927 + 2.422723,, _i — l .921 142„_2+0. 47949 2„_ 3 (7) 


This equation may be regarded as a recurrence formula; the 
value z„ may be interpreted as the deflection to come and 
may be found easily from the three preceding deflections 
Zn-u 2 , 1 - 2 , and z n _ 3 - Then with the newly found value z n 
and with 2 „_, and 2 „_ 2 , the next deflection can be found, 
and so on. This process thus gives a step-by-step deriva- 
tion of the time history of deflection and may be carried 
out as far as is desired. Of course the process must start 
with known initial values of z. These values can be found 
with the aid of the initial conditions of the problem by 
means of the following approach. 

The difference equations for the first and second deriva- 
tives at the third ordinate of four successive ordinates are 
(see appendix B) 

w„ = g^-(2w„ +1 +3w„ — 6w„_i+in„_ 2 ) 

w n =^r i {w n+ i—2w n -\-w n - l ) 

If these equations are taken to apply at t ~ 0 {n— 0), they 
become 

m,=-^(2w 1 + 3wo-6®-i+W-2) (8) 

w 0 =- ~ (W! — 2 w 0 + i) 


(9) 
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For the problem under consideration the primary initial 
conditions are that, at t= 0, the displacement and velocity 
are zero. By use of equation (1) or by reasoning from 
Newton’s second law, a secondary initial condition can be 
established — that is, the acceleration immediately following 
the application of the unit force must be 1/m. In equation 
form these conditions are 

w 0 =0 

w 0 =0 

1 

w Q =— 

m 

By substitution of these values into equations (8) and (9) 

w 

and by use of the notation z= -p]k’ following relations 
can be found to exist between the ordinates: 

2 0 =0 'j 

3_2=0. 24-82! i- (10) 

s_! = 0.04 — z i J 

Substitution of these values into equation (7), with n set 
equal to 1, gives an equation from which z u the deflection 
at t=e, may be evaluated. Three successive deflections 
can now be established: the deflection at £=e, the zero 


step-by-step evaluation of succeeding deflections proceeds in 
a straight-forward manner — that is, equation (7) is now 
evaluated for n= 2, then for n— 3, and so on. The response 
of the oscillator for the physical constants listed previously 
is given in figure 2. The comparison between the difference 
solution shown in this figure and the exact solution of 
equation (1) is seen to be good. As a matter of interest, 
the solution is also shown in this figure that is obtained by 
the use of the parabolic end-ordinate derivative which 
involves only three successive ordinates. The agreement in 
this case is seen to be quite bad. If equations (2) and (3) 
had been used, on the other hand, the difference solution 
(in this case for w„ + 1 ) would correspond to that given for 
the cubic end-ordinate derivative. 

RECURRENCE MATRIX EQUATION FOR RESPONSE OF AN AIRPLANE 

IN A GUST 

In order to help the reader to obtain a perspective of what 
is to be covered in this section, the following basic phases of 
the analysis are given: 

(1) The assumptions are stated. 

(2) The coordinate system and displacements are defined. 

(3) The aerodynamic lift and moment are defined. 

(4) The normal and torsional dynamic loadings (inertia 
forces, aerodynamic force s, and fuselage forces) on the wing 
are derived. 

(5) The equations of elastic deformation — wing vertical 


deflection at t= 0, and a fictitious deflection for t—~e given 
by equation (10). In the real problem no deflection exists 
for t less than zero; the assumption that a deflection does 
exist before t is zero is simply a means for allowing the recur- 
rence formula, equation (7), to apply at the origin as well 
as at later values of time. Furthermore, no violation is 
made of the conditions under consideration because, mathe- 
matically, the response after f=0 is not influenced by the 
response that may exist before <=0, so long as the initial 
conditions are satisfied. The process just described for 
treating the initial conditions is actually not different from 
the procedure commonly employed in difference-equation 
approaches, in which exterior points near a region under 
consideration are written in terms of the interior points by 
means of the boundary conditions. 

With the initial values of deflection thus established the 


motion, wing rotation, and fuselage bending — are given. 

(6) The dynamic loadings on the wing are transformed 
into difference equations. 

(7) The equations of elastic deformation and the difference 
equations for loading are combined to give the recurrence 
matrix equation for response. 

In succeeding sections the initial response is determined, 
the method for evaluating the gust forces is shown, and the 
method for computing the loads and stresses is indicated. 

Assumptions. — In this analysis an attempt is made to 
obtain the simplest and most direct solution to the problem 
with a minimum of simplifying assumptions. The case 
treated is that of an airplane having an essentially straight 
wing and penetrating a gust of known structure. The tail 
is considered to penetrate subsequently the same gust as 
does the wing. The assumptions made are as follows: 



Assumptions pertaining to elasticity and airplane motion: 

(1) The usual assumptions of engineering beam theory are 
made. 

(2) The fuselage is free to pitch and move vertically. The 
portion of the fuselage in front of the elastic axis of the wing 
is assumed for convenience to be rigid. The portion of the 
fuselage rearward of the elastic axis is assumed flexible, and 
the elastic deflection is assumed to be given by a single 
modal function. 

(3) The tail is assumed rigid. 

Assumptions pertaining to aerodynamic forces: 

(1) Aerodynamic strip theory applies. Three-dimensional 
effects, however, may be taken into account approximately 
by means of over-all corrections. Some such corrections are 
indicated. 

(2) The gust force and forces due to vertical and pitching 
motion are the only tail forces considered. Other forces of 
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(3) Aerodynamic forces on the fuselage are neglected. 

Coordinate system and displacements. — Position on the 
airplane is denoted by an orthogonal system of axes.. The 
origin is at the intersection of the wing elastic axis with the 
plane of symmetry of the airplane: the w-axis tuds positive 
upward, the z-axis runs along the fuselage positive in the 
rearward direction, and the y-axis runs spanwise. The wing 
semispan is considered to be divided into six, not necessarily 
equal, sections, with a station point at the middle of each 
section. (See fig. 3.) More or fewer stations could be 
chosen, but it is believed that six is a fair compromise be- 
tween the amount of labor involved in setting up a solution 
and the accuracy desired. The interval between stations is 
designated by the number of the station at the outboard end 
of the interval. Station 0 is located near the wing root and 
generally may be located where the fuselage intersects the 
wing. In this way the concentrated forces of the fuselage 
are allowed to act through station 0. The other five stations 
are then located in any convenient manner so as to fall at 
concentrated mass locations or at points which represent the 
average of distributed masses, station 5 being nearest the tip. 
The total mass within a section is assumed to be concen- 
trated at the station point, and the average of the section 
geometry (chord, elastic axis position, and so on) is assumed 
to apply. In this way the wing is assumed to be a beam 
subject to six load concentrations and as such will have a 
linear moment variation between each station. The further 


assumption is made that the 


1 

El 


variation is linear between 


each station. With these assumptions for the El variation 
and concentrated load locations, equations for deflection at 
each station point may be derived (appendix C) by direct 
analytical treatment. 

The displacements of the cross section at each station of 
the wing are given as the deflection of and rotation about the 
wing elastic axis as shown in figure 4. The fuselage dis- 
placements are shown in figure 5 and are given by the 
equations: 


Wf=w{ 0) — <p(0)a; 


( 11 ) 


for the forward section and 


w f =w(Q)~- <p{0)x + Wyb (12) 

for the rearward section. The function IT, is taken as the 
fundamental mode of vibration of the fuselage and tail 
assembly, when the fuselage is considered to be clamped as a 
cantilever beam at the elastic-axis location of the wing, and 

is given in terms of a unit deflection at the ^ -chord position 

on the tail. With this function to represent the elastic 
deformation of the fuselage the deflection and angle of 
attack of the tail are found with the aid of equation (12) to be 


at 


_ dw f ~ | 


where 


w(0) — ^(0)x, + 5 

(13) 

-1 =<KO)-50! 

(14) 

Jx=x t 

_dwr\ 
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Figure 3. — Division of wing into sections. 




Aerodynamic lift and moment. — Before going into the 
details of the analysis it is felt worthwhile to define and 
describe the nature of the aerodynamic forces to which the 
wing is subjected. These forces originate from two sources: 
they arise directly from the gust encountered, and they arise 
from the ensuing airplane motion. The equations for the 
aerodynamic lift and moment that develop are herein set up 
in a convenient form on the basis of work given in references 
1 to 4. In these investigations various methods for separat- 
ing the lift forces have been used, but the particular method 
for separating these forces is not important so long as they 
are taken into account properly. 

In the present report the aerodynamic lift and moment 
are considered to be composed of two parts: one part, desig- 
nated as the lift or moment due to circulation, which includes 
all lift forces or moments exclusive of aerodynamic inertia 
effects and the other part, which is due solely to these inertia 
effects. These lift forces and moments can be resolved into 
the force systems acting on the airfoil as shown in sketch 1 
(forces due to circulation) and sketch 2 (inertia force and 
moment) . 
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The force L g is the lift force developed by the gust and 
corresponds to the gust force that would develop on the 
airfoil considered rigid and restrained against vertical motion. 
All the other forces occur as a result of motion of the airfoil. 
These forces, as well as the gust force, are given for an 
interval l of the span by the equations: For the forces due 
to circulation, 

Lg—m A irpclU P IP ^(t— T)dr (15) 

Jo or 

Li= m A irpclU ^ [”[/£ — w-\-c [1 — (f~ r)] dr 

(16) 

Z 2 =— u<p (17) 

and for the inertia force and moment, 

L s =^~ (18) 

= (19) 

where 

m A factor which can bo used to introduce over-all com- 
pressibility and aspect-ratio corrections; in this 

report the factor is assumed to bo given by 

A 

2 + A^/I—M 2 

1 — F lift function which denotes the growth of lift on 
an airfoil following a sudden change in angle of 
attack 

\p lift function which denotes the growth of lift on a 

rigid airfoil entering a sharp-edge gust 
The functions 1 — F and p and the correction m A are 
established as follows. In reference 5, approximate equa- 
tions are derived which give the lift-coefficient form of the 
growth of lift on a finite wing following a sudden change in 
angle of attack or due to the penetration of a sharp-edge 
gust. The equations may conveniently be considered as the 
product of a factor, which may be regarded as a lift-curve 
slope, and an unsteady lift function, designated by 1— $ for 
the function due to the angle-of-attack change and by \p 
for the function due to the sharp-edge gust. These unsteady 
lift functions are shown in figures 6 and 7 and are given by 
the following equations: For the 1— $ functions 



(1 — F)^3= 1 — 0.283e-°- 540s 

(20a) 


(1— F) x _ 6 =l — 0.361e~°- 381s 

(20b) 


(1— F) x= „ = ! — 0.1 65e -0 045s —0.335e _0 - 300s 

(20c) 

and for the p functions 



^ =3 = 1 -0.679e-°- 558s -0.227e- 3 - 20s 

(21a) 

V Pa 

„ 6 = 1 — 0.448e~°- 290s — 0.272e~°- 726s — 0.193e _3 00s 

(21b) 

Pa, 

: „ = 1— 0.236e-° 068s — 0.513e-°- 364s — 0.171e- 2 - 42s 

(21c) 



Figure 6. — Lift functions for sudden change in angle of attack. (Sec equations (20).) 



Figure 7.— Lift functions for wings entering a sharp-edge gust. (Sec equations (21) 

and (22).) 

Equations (21) are based on equations of reference 5 ; whereas 
equation (22) is the p function that is suggested for wings of 
infinite aspect ratio in reference 3. Inspection of equations 
(20) shows that the $ function for aspect ratios 3 and 6 is 
given by a single exponential term. It is probable that the 
$ function for higher aspect ratios, say 10 and even 20, 
may also be given to a sufficient approximation by a single 
exponential term. Therefore, the assumption is made that 
in general F may be represented by an equation of the form 

F=a,<r* s (23) 

Interpolation, for example, of the curves in figure 6 shows 
that the F function for aspect ratio 10 might be approximated 
by the equation: 

F A „ IO =0.41e-°- 3s (24) 

The analysis does not necessarily limit F to a single ex- 
ponential term. Other terms could be added with some 
increase in labor, but the degree of refinement obtained 
is not expected to add much to the over-all accuracy of the 
solution. 

Although the functions given by equations (20) to (22) 
are known to approximate the true functions quite well over 


7 


A RECURRENCE MATRIX SOLUTION FOR THE DYNAMIC RESPONSE OF AIRCRAFT IN GUSTS 


equations (21) do not vanish, as they should, when i=0. When 
used in the computational procedures given hereinafter, these 
functions, therefore, are to be taken as zero when t= 0. 
Another point to note is that the variable s is given in terms 
of a reference chord c 0 ; thus this variable as applied to the 
wing is different, in general, from the variable as applied 
to the tail. 

Examination of the values of lift-curve slope, which were 
stated to be present in the equations taken from reference 5, 
reveals that they may be approximated with good accuracy 
by the product of 2ir and the often-used aspect-ratio cor- 

A 

rection f° r sccac b ; incompressible flow. In the present 

report it is assumed that compressibility and aspect-ratio 
corrections can be made by replacing this aspect-ratio 
correction by a compressible aspect-ratio correction defined 

by 3 


r + 2 


where A'—A-^l—M 2 , and by multiplying this 


correction by the Glauert-Prandtl Mach number correction 

, 1 to give the product m A . The procedure then for 

VI — M l 

taking into account three-dimensional and compressibility 
effects in the present analysis is to determine m A from the 
forward speed and aspect ratio of the wing and to use the 
1 — <f> and \p functions, equations (20) to (24), for the aspect 
ratio which is nearest that of the wing. 

Some word of explanation of equation (16) might be 
worthwhile at this point. The <f>(< — r) function is associated 
with the lift forces which are due to the wake. Without this 
term the equation would yield the steady lift corresponding 
to the instantaneous values of angle of attack and vertical 
velocity. If equation (16) is integrated by parts and the 
conditions are stipulated that w, w, <p, and ip are all zero at 
t = 0, the following equation may be found: 


L x = f}cl<b 0 w—(\ — ® 0 )pclw+pclU^l 

(1 — <f> 0 )/3c 2 Z^— + w $(t— T )d T — 

PclU J y>4>(f— t)(It — !)/: ^(f-r)dr (25) 

where /3 has been introduced as a forward-speed and aspect- 
ratio parameter defined by the equation 

f} = m A TpU (26) 

With reference to equation (23), 4> 0 and 4> 0 in equation (25) 
would have the values 


i 2U . 

<f> 0 = — — 

Cq 

The form of L x given by equation (25) is presented because 
this form is more convenient to use in the present analysis. 

For this analysis the total lift and moment acting at the 
elastic-axis location are desired. For the present, the total 
lift L and moment M of the forces due to circulation are 
found; the inertia force and moment are to be treated 


separately. Summation of all the lift forces due to circula- 
tion and summation of the moments of these forces about 
the elastic axis gives the desired equations for the aero- 
dynamic lift and moment acting on the airfoil over an 
interval l as follows: 

i=il + i 2 + ig (27) 

M =( a “£) il- (T -a)z 2 +(a-£)4 (28) 

The loading on the wing. — The normal and torsional 
dynamic loads that are held in equilibrium by the elastic 
restoring forces of the wing may be found by considering all 
the aerodynamic and inertia forces that act on the wing. 
The mass situated at any station (see fig. 4) can be shown to 
have an inertia normal force equal to 

— mw-\-melp 

and an inertia torsional moment about the elastic axis equal 
to 

meiv—mk 2 if 5 

If the aerodynamic forces and moments (see equations (18), 
(19), (27), and (28)) are added to these inertia loadings, the 
total normal and torsional loadings on the wing at each 
station are found to be given, respectively, by the equations: 

p= — mw+m«p+I+ L 3 

q — mew — mk 2 ip-\-M—(^—a^ L 3 -\-M x 

The terms L 3 and M x ordinarily would appear with the aero- 
dynamic lift and moment values but are treated separately 
so that they can be combined with the structural mass 
terms. If use is made of equations (18) and (19), the load- 
ing equations become 

p= — mw + rne<p-\-L (29) 

q = meib—mk 2 (p-\-M (30) 

where 

— ( i Trplc 2 \ 

m=( m -\ — - — 1 

«+nr(|-f)] 

The terms appearing with the structural mass quantities in 
the definitions of m, me, and mk 2 are the terms which are 
commonly associated with apparent mass effects. 

The value of the shear forces L f and the moment M f trans- 
mitted to the wing by the fuselage can be found in the 
following manner. From equations (11) and (12) the values 
of the inertia loading on the forward and rearward sections 
of the fuselage can be shown to be given, respectively, by 
the equations: 


p /= — m , [ia(0) — «5(0)x] 

(31) 

p f = — m f [ii>(0 )— v5(0)x + '5IFi] 

(32) 
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Integration of these inertia loadings over the length of the 
fuselage and addition of the aerodynamic tail load 2 L t give 
the value of the total load transmitted to the wing; one-half 
of this load is designated by L f and is assumed to act at 
station 0, the other half being considered to act through the 
corresponding station on the other half of the wing. Inte- 
gration of the moment of the inertia loading about the 
elastic-axis location and addition of the moment —2x t L, of 
the tail forces give the total moment due to the fuselage; 
one-half of the moment is designated M f and acts at sta- 
tion 0. The values of L, and M, thus found can be given 
by the equations: 

Lf— — Miw(0) +M 2 v?(0 ) - — M 3 5 + L, (33) 

M f =M 2 w(0)—M i 'ip(0)+M 3 d — x t L t (34) 

where the M/s are considered to be generalized masses 
defined as follows : 


form, the bending and rotational displacements are related 
to the normal and torque loadings by the well-known 
expressions : 


b 2 w_ 

by* EI by 2 P 


(36) 


by by 


= 2 


(37) 


where in this instance p and 5 are the loadings per unit 
length of beam. In addition to these two equations which 
are considered to apply to the wing, an equation for com- 
puting the elastic deformations of the fuselage may be found; 
this equation may be found in the following manner. The 
rearward part of the fuselage is considered to be a cantilever 
beam subjected to the inertia loading given by equation 
(32) and the tail force 2 L t . If equation (36) is applied 
to the fuselage and use is made of equations (12.) and (32), 
the following equation for fuselage bending results: 


Mi 2 ^ 

f X ‘m f dx 

/ -t-n 

M 2 =i 

I rrifxdx 

4 

^ i <N 

II 

rx t 

| W/lTi dx 

!o 

fie*' 

II 

£ 

C x ‘ 

j m f x 2 dx 


j" dx 

M e =i 

f m f W 2 dx 

Jo 


The generalized mass constant M s , although not appearing 
in equations (33) or (34), is included in this group because 
it occurs in a subsequent part of the analysis. In the 
derivation of equation (34), the aerodynamic moment of 

the tail about the tail --chord position is neglected since 

it is considered to be small in comparison with the value 
x t L t . The lift on the tail L t can be found by application 
of equation (27) to the tail surface. In this case the >t> 
function appropriate to the tail should be chosen and the 
values of displacement w and p should be replaced by 
Wf{x t ) and a t , the values of deflection and angle of attack 

at the tail i-chord position. These values are given by 

equations (13) and (14). 

Matrix equation of equilibrium. — The problem of comput- 
ing the response may be considered to be one of the deter- 
mination of the deflection and rotation of a beam which is 
subjected to normal and torque loadings. In differential 


ft 2 a2tt/ 

t^CO) — v5(0)x + S ^ + 2 (38) 


in which L t must be treated properly as a concentrated load 
and If is the bending moment of inertia of the fuselage. 
Since W, represents a vibration modal function, the following 
relation exists: 


b' 2 2T „ 


where w, is the frequency of vibration associated with IT, . 
Equation (38) may therefore be written 

8mfa>/Wi = — 7 n r [w(0) — vi(0)r + 5 W T i] +2 L t 


Multiplication of this equation through by IT; and integration 
between 0 and x t results in the following equation for fuselage 
bending 

Uf 2 M$d= — J\d 3 ib(0)-}-M 3 'ip{0) — M§8 T Lt (39) 

where M 3 , M s , and M 6 are defined by equations (35). 

Equations (36), (37), and (39), when the loadings given by' 
equations (29) and (30) are consideied, are seen to be rather 
involved integrodifferential equations but describe com- 
pletely the motion of the airplane. The problem is to find 
functions w, <p, and 8 which satisfy these equations and which 
satisfy both the boundary conditions and the initial conditions. 

The problem of finding the w and <p functions may 7 be 
simplified considerably by reducing the rather complicated 
equations of motion to a simplified and systematic algebraic 
form. The first step (see appendix C) is to replace the dif- 
ferential equations (36) and (37) for wing deflection and wing 
rotation by the following simple matrix equations: 

[A] \w\ = \j)\ (40) 

[B] M = kl (41) 
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The matrices in these equations are defined in appendix C 
(see equations (C22) and (C23) and equations (C29) and 
(C30), respectively) and have been derived on the basis that 
the displacements along the semispan are given at six stations. 

Equations (40) and (41) and the fuselage deflection coeffi- 
cient d are now combined in a single matrix equation of the 
form indicated as follows: 


r~ 

© 

o 

o 

l_ 


5 


0 

0 [A] 0 


M 

= 

bl 

i 

a 

o 

o 


1 <e\ 


Iffl 


This form is chosen because it will be useful subsequently. 
With the notation given in appendix A, equation (42) may be 
abbreviated to the form: 


[C] 


= \P\ 


(43) 


This equation may be regarded as the loading matrix equa- 
tion of equilibrium; it relates the loadings to the displace- 
ments by linear simultaneous equations. The boundary con- 
ditions are automatically satisfied when this equation is used 
because they had to be taken into account when the sub- 
matrices [A] and \B\ were derived. Only the initial conditions 
remain to be satisfied and these are treated separately in a 
subsequent section. 

Transformation of the loading equations into difference 
form.— The loading equations are now simplified by replacing 
the time derivatives by difference equations. If equation 
(5) is used to replace the derivative in equations (29) and 
(30), the values of the loading at the nth time interval are 
found to be 


Pn = —f (2w»-5w„. 1 + 4w»_ 2 -w„_ 3 ) + -^ (2<p n — 5y„_i + 
^‘Pn-2~<Pn~i)-\-L n (44) 

7 YIC 771 Jc ^ 

q.n=—£ (2w,-5w»_ 1 +4w„_ 2 -w„_ 3 ) -j- (2<p n — 5<p„_i-i- 

4<p K _ 2 — <Pn-3) + M„ (45) 

The values L„ and M n are found by determining the expres- 
sions for Li, Li, and L e at t=nt (see equations (27) and (28)). 
Of these L y is the most complicated, since it (see equation 
(25)) involves three unsteady lift integrals of the Duhamel 
type. Fortunately, however, a rather simple recurrence 
relation can be developed which allows the calculation of the 
value of these integrals at a given time interval directly from 
the value at the previous time interval. This derivation is 
presented in appendix D and is made possible because the 
<l> function is of an exponential form. (Where the 4> fimetion 
is given by more than one exponential term, a recurrence 


relation for each term may be written.) From the derivation 
in appendix D, therefore, the value of the three integrals at 
the Tith time interval may be given as follows: 

F n +~- $0 [W« (f — c) So] (4 6) 

where 

F n = e ~f‘ F n _ i +gw n _ i + g‘ '<p n _ i 

in which g and g’ are defined by equations (A5) in appendix 
A. With this expression to replace the value of the integrals 

in equation (25), the value of L, n may be written 
% 

L ln = f3cl w„ — (1 — 4> 0 ) ficlw n + 

Pel [t/(i -*„)-^$ 0 -(<i. 0 +f $o)c “)] <p n + 

(1— f>„)/3c 2 z(!^)^ + F„ (47) 

With the use of difference equation (4), this equation may 
be transformed finally into the form: 

L Xn = d 0 w„ + diW n _ i + d 2 w n _ 2 + d 3 w„ _ 3 + d 0 '<p„ + 

di <p n -i-\-di <p n -i-\-di ipn-s~\~F n (48) 

where the d’s are defined in appendix A. Likewise, from 
equations (4), (17), and (26), A 2 „ may be written 

(lllpn— 18^„_ I 4-9(p„_2— 2<p„_ 3 ) 

n Z4€ 

If L~ , L, , and the value L„ of the gust force at t = nt 
are introduced into equation (44), the value of p at the nth 
time interval can be shown to be given by the equation: 

P n = yoWn + yiWn-\+n2W't-2+ri3Wn-3+yo'<P l L-i- 

V\ ‘Pn-\ J rV 2 <Pn- 2 J r V3 '<Pn - 3 + F n + Z g(; (49) 

where the rj’s are coefficients which arc given by equations 
(A3) in appendix A. In a similar manner, the equation for 
g (equation (45)) can be reduced to the form: 

q n = v 0 w n + Vi W n _ 1 + v 2 W n _ 2 + v s W n _ 3 + roV* + "1 V* - 1 + 

r 2 Vn-2+J , 3 Vii-3+(^« — 4 ^ F n -\-^a— L gn (50) 

where the r’s are given by equations (A4) in appendix A. 

The value of aerodynamic lift acting at the tail “-chord L t 

is found most conveniently by applying equation (49) to 
onc-half of the tail surface. This application is made by 
modifying the y coefficients as follows: The mass value m 

is set equal to zero, — is taken as c is replaced by c ,, and 
C 4 : 
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pci is replaced by P„ defined as the forward-speed and 
aspect-ratio parameter of the tail by the equation: 

m Al T P US, (51) 

In addition, w and <p are replaced by the deflection and 
rotation of the tail given by equations (13) and (14). With 
these substitutions the value of L, is found to be 

r* n 

4„ =f 0 w(0 ) n +/iw(0)„_ 1 +/2'w(0)„_ 2 +/ 3 'U>(0), i _ 3 +/oV(0), I + 

/lV(0)n-l+/2V(0)n-2+/3V(0)7!-3+/o5 re +/l5 B _i-l- 

i + 4, (52) 

n 

where 

F 1 t n =^~ y,e Ft n _ 1 -\-jw(0) n -i-\-j , <p(0) n -i+j8 n -i (53) 

Lg t is one-half the tail gust force at t=nt and the/’s and j’s 
are defined by equations (A7) and (All), respectively, in 
appendix A. 

With equation (52) and difference equation (5), equations 
(33) and (34) for L f and M f and equation (39) for fuselage 
bending may be reduced readily to the following form: 

A/„ = 7oW(0) re +7iW(0) rc _i-|-Y2W ( 0 ) 71-2 + 73 ^ (0) n -3 + 7oV (0)n + 

7lV(0)n-l+72V(0)n-2 + 73V(0))l-3 + 7o5j! + 7 i5ti-1 + 

725,7-2 + 735,7-3 + 4,,+ Lg t ^ (54) 

Mf n — HoW (0)„. + Hi'!e(0) n _ l + jU2'W (0)m-2 + / i 3'W (0) K _3 + /io , ^’(0)„ + 
MlV(0)n-l + M2 / ¥ 3 (0) n -2 + M3 , ¥ 5 (0)7!-3+Mo5, l + Ml5 n -l + 

MaS, 7-2 + 1*35,7-3 t n — x t Lg t ^ (55) 

0=r 0 w(0) re +r 1 w(0) B _i+r 2 w(0)„_2+r 3 w(0) B -3+r o V(0) B + 

/’iV(0)„_i + r 2 V(0) n -2 + ?'3V(0) n _3 + / , o57t + /‘i5„_i + 

r 2 8 n -2 + r 3 Sn -3 + F ln + L gtn (5 6) 

where the 7’s, +s, and r’s are given by equations (A8) to (AlO) 
in appendix A. 

The complete set of loading equations as well as the 
fuselage bending equations are now available in difference 
form. Equations (49) and (50) apply at each spanwise 
station and in addition the values of L f and M r must be 
introduced at station 0. The coefficients q, v, 7, and so 
forth are seen to involve only the physical properties of the 
airplane structure, the forward-speed and aspect-ratio 
parameters given by equations (26) and (51), certain con- 
stants derived from the unsteady lift function, and the time 
interval. The time interval e that is chosen should be 
fairly small in comparison with the natural period of the 
fundamental mode in bending of the wing. To serve as a 
guide an interval chosen near 1/30 of the estimated period 
of vibration of the fundamental mode appears to be quite 
satisfactory. Of course, some caution should be observed 


in the choice of this interval if the airplane is near a critical 
condition where some mode other than the fundamental 
may predominate. For example, if the airplane is flying 
near the flutter speed, the characteristic frequency of the 
response may be near the natural torsional frequency of the 
wing. The time interval should be modified accordingly. 

Recurrence matrix equation for response. — Equations (49), 
(50), (54), and (55) for loading, equation (56) for fuselage 
bending, and the equilibrium equation (43) may now be 
combined to give the recurrence matrix equation for re- 
sponse. In order to simplify the process of combining these 
equations, only the abbreviated or symbolic form of the 
matrices which occur are used. The definitions of these 
matrices are given, unless otherwise stated, in a complete 
group in appendix A. 

Application of equations (49) and (50) to each of the span- 
wise stations and of equations (54) and (55) to station 0 leads 
to a set of loading equations which may be put in the matrix 
form given by the following equations: 


I p\n — |7ol 5 n + 1 7l 1 5 t, _ 1 + I 72 I 5„ _ 2 + 1 73 ! 5 ra — 3 + [ 770 ] I'W | ra + 
[)ji]|w|„-i + [r;2]|wU-2+[’?3]|'^U-3 + [’?o / ]|^|77 + 

h.lkl.-i+MM.-.+MI#»U-3+|l^l+IAl| + 
111 (4+4,)* 

[ S Ira == 1 Mol 5,1 + I At 1 ! 5,, — 1 + I /x 2 | 5 t ,_ 2 + I M 3 I 5t, _ 3 + [>>o] 1^1 ra + 

[''l]|w|„-i + [r 2 ] 1^177-2 + [^\w\n-a-\-[vo]\<p\n + 

Mt 7-1 + [»'2 , ]|¥>|77-2+ h3']|#>U-3 + £& — |4| + 

+ \x t \(F t +L gt ) n 


L"1 J 

141 


where 


\F\ n =e 7< |E| B -i + [g]|w|, ! - 1 + [g , ]!<p|,,-i 


(57) 


(58) 

(59) 


(E,)„ = e y ‘ e (E,) n - 1 +jw(0) n - 1 + i )V (0 )ti_i + jS n -i (60) 


Equations (57) and (58) and equation (56) may now be 
combined to form the following matrix equation: 


0 


~ 4 

ho] 

[4]" 


5 


\v\ 

= 

Ito| 

ho] 

ho'] 


M 

+ 

hi 

n 

Jmo 

ho] 

ho'L 


\<p 1 

n 


4 h 2 J h 2 'J 
I72I ha] ha'] 
LI/I2I [r 2 ] [4]J 
1 0 

ll [/] 

[“-I] J 


r 1 In] In'] 
7.1 hi] h/] 
Tii I h>] [+]j 
"4 [r 3 \ In'] 

1 73 1 ha] ha'] 
n- 2 L|4| [ ^3] [ ^3'] J 

4 + 4, 


+ 


14+141 1 


+ 


+ 


( 61 ) 


A RECURRENCE MATRIX SOLUTION FOR THE DYNAMIC RESPONSE OF AIRCRAFT IN GUSTS 


11 


For simplicity, this equation may be abbreviated to the 
form: 



8 


s 


8 


8 

\P\,-[S 0 ] 

w 

+ [&] 

w 

+[Sd 

w 

+ [S 3 ] 

w 


<P 

re 

<f> 

71—1 


re — 2 

<P 


rwi+iz, 


where 


\F\ n =[e]\F\ n ^+\Wi 


s 

w 

<P n - 1 


(62) 


(63) 


and the matrix \L g \ n is defined in the section entitled “Der- 
ivation of Gust Forces.” 

Substitution of equation (62) in equation (43) gives 



8 


8 


5 


8 


8 

[C] 

w 

= [S 0 ] 

w 

+ [£,] 

w 

+[s 2 ] ! 

w 

+ [S a ] 

w 


<p 

n 

<P 

re 

<P : 

re — 1 

•f 

re — 2 

<P 


+ 


r^J 


IM+I l s 


(64) 


With the use of the notations 


[D]=[[G]-[S„]] (65) 

and 




8 


8 


8 




1 Q\.= 

[S 1 ] 

w 

+ [SJ 

w 

+ [S 3 ] 

w 

+r^J 





<P 

re— 1 

<P 

re —2 

<P 

re —3 


re 


equation (64) may be written simply 


( 66 ) 


5 


From equation (68) the complete response of the airplane 
can be computed once the character of the gust is known. 
The matrix of gust-force values |Z*| n can be determined by 
the procedure given in the section entitled “Derivation of 
Gust Forces.” The initial conditions (treated in the follow- 
ing section) are used to obtain three successive initial sets of 
the displacements. With these sets of displacements the 
next set may be obtained by application of equation (68). 
With the newly found set and the preceding sets of displace- 
ments, the next set may then be found, and so forth, until a 
sufficient time history of the displacements is found from 
which maximum loading conditions may be determined. 

The reason for using the difference form of the derivatives 
as given by equations (4) and (5) might now be given. 
Equation (64) may be considered a differential equation, 
since the matrix [C] contains the spanwise derivative matrices 
[.4] and [B] and may be likened to the differential equation 
which relates the load to the deflection for a beam. The 
unknowns are the deflections at time n. The right-hand 
terms correspond to the loading, the first term being the only 
one which is not known since it contains the unknown de- 
flection. The subsequent inversion of the matrix [D] leads 
to, in effect, the solution to this differential equation and, in 
the beam analogy, corresponds to the integration of the 
loading to obtain the deflection. When numerical methods 
are used, the deflection may be computed with good accuracy 
by integration of the loading. On the other hand, if the 
difference equations which apply at an interior ordinate 
had been used, the matrix [ C ] would have appeared as an 
operator on one of the known deflections on the right-hand 
side of the equation. Effectively, its operation would be to 
differentiate a known deflection,. and in the beam analogy 
this operation corresponds to the attempt to obtain the load 
which caused a given deflection. This process, however, 
cannot be done with accuracy when numerical methods are 
used because of the difficulty encountered in the form of 
small differences of large numbers. The difference equations 
which apply at an outer ordinate and, consequently, lead 
to an integration process, therefore, have to be used. 


[D] 


w 


= \Q\n 


<P\n 


(67) 


Multiplying through by the reciprocal of [ D ] gives finally 
the equation 

i SI 


= [D]~'\Q\ n 


( 68 ) 


This equation gives the displacements that apply at time n 
in terms of the displacements that occurred at several pre- 
ceding values of time (see equations (63) and (66) for the 
definitions of ]p| M and |Q|„. 


DERIVATION OF THE INITIAL RESPONSE 

As has been mentioned, some initial values of deflection 
are needed before equation (68) can be used. This section 
shows how these values are obtained. The airplane, just 
before gust penetration, is considered to be in level flight, 
and all displacements which follow are given relative to this 
level-flight condition. The initial conditions are that the 
vertical displacements, vertical velocity, wing rotation, and 
angular velocity are all zero. The gust force can be shown to 
start from zero and, therefore, by Newton’s second law the 
additional initial condition can be established that the 
acceleration must be zero at the start of the response. These 
conditions can be satisfied, and the beginning of the response 
can be found by means of the analysis which follows. 
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The initial response is assumed to be given in terms of 
four successive ordinates, denoted by W- 2 ,W- i,tv 0 , and Wijthe 
w 0 ordinate is considered, as in the case of the damped oscil- 
lator, to locate the origin of time. The first and second 
derivatives at the w 0 ordinate are given by equations (8) 
and (9), respectively. By virtue of the initial conditions 
(the vanishing of the deflection, velocity, and accelerations 
at t= 0), the ordinate w a and the derivatives given by equa- 
tions (8) and (9) must be zero; therefore, the ordinates w_ 2 
and w—i are found to be related to the ordinate w x by the follow- 
ing relations : 

w_2=— 8Wi (69) 

w_i = -Wi (70) 


These relations are general and must apply for deflection 
and rotation at each of the spanwise stations and for the 
fuselage deflection as well; that is, the displacements at 
t— — 2t must be minus eight times the displacements at t=e, 
and the displacements at t == — e must be the negative of those 
at t=t. Substituting these conditions in equation (64), 
taking n as equal to 1, and using the condition that the dis- 
placements are zero at t— 0 give the following matrix 
equation in terms of the displacement at t—e alone: 


[[£>]+ [«y + s [s 3 }] 


8 

w | 

<P 


=m\L,i 


(71) 


The term j/ <T | 1 is zero and therefore does not appear in this 
equation. Solution of this equation gives the values of the dis- 
placements that occur at t=e (n= 1 ). 

The three sets of initial displacements required to proceed 
with equation (68) are thus known: the set of deflections 
found att=e, the zero set at t=Q, and the set at t= — e given 
by equation (70), or simply the negative of the displacements 
which were found at t—t. In the actual case no displace- 
ments are present at t= — e, but these displacements may be 
regarded as being of a fictitious nature the only purpose of 
which is to allow the step-by-step evaluation of the dis- 
placements to be started easily. 

DERIVATION OF GUST FORCES 

The matrix \L g \ n which appears in the response equation 
(68) is derived as follows. From equation (15) and the 
notation of equation (26), the total gust force acting over a 
station section at the nth time interval may be given by the 
equation 

L gn =pcl £ t(ne-T)dr (72) 


with this response the response to any gust may be found 
directly by superposition. 

Thus for the case of a sharp-edge gust, equation (72) 
reduces simply to 


L Sn —&clv 


(73) 


where \p n is the value of the 4/ function at t—nt. This 
equation when applied to each of the spanwise stations leads 
directly to the matrix equation for gust force: 


\L g \ n — P 


ColoVf) 

CihVi 

C 2 l 2 V 2 

C3I3V3 

C4I4V4 

C5I5VS 


<Pn 


(74) 


If the gust is uniform in the spanwise direction, the w’s in 
this equation will all be equal. 

In a similar manner, one-half the gust force acting on the 
tail due to a sharp-edge gust may be shown to be 

I's, n = PtVo'l't n (75) 

where the gust gradient is assumed to be the same as for 
station 0 and 4 / t n is the value of the function for the tail. 
This equation and equation (74) may now be combined to 
give the desired matrix \Lj\ n as follows: 


(76) 


In the application of this equation it should be kept in mind 
that L et does not begin to act until the tail starts to penetrate 
the gust. The time interval at which penetration starts 

may be taken as the interval nearest to the quantity 

COMPUTATION OF LOADS AND STRESSES 

Once the time history of the displacements has been found 
from equation (68), the normal or torque loading on the 
wing can be found with little additional work. If the nota- 
tion of equation (66) is used, equation (62) may be written 


1 _ 

L' t 

rt 

\L e \ 


hV 0 

0 

0 

0 

0 

0 

0 


0 

PcJoV 0 
pC\hv x 
Pc 2 l 2 V 2 
PC 3 l 3 V 3 
PC^liVi 
Pc 5 l 5 V s _ 


The integral in this expression is also of the Duhamel type 
and since the ^ function is expressed by exponential terms 
(see equations (21) ), the integral may be evaluated quickly 
by a method similar to that developed in appendix D. The 
procedure of computing the gust force by this equation 
and then the response is not recommended, however, since 
a complete response evaluation would have to be made for 
each yust considered. Instead the procedure 'eeom xercer 


\P\n = [S,] 


+ I Q\n 


(77) 


This equation shows that the loading matrix \P\ may be found 
by adding an easily computed matrix to the matrix |Q|, the 
value of which will have already been determined in the 
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to be defined in terms of the normal and torque loadings, and 
either of these loadings may be found independently of the 
other. 

The loadings thus found are considered to be applied 
statically, and the stresses are then found by ordinary static 
means. Since the loadings can be computed for any value 
of time, the complete stress history of any point in the 
structure may be computed. In general, the maximum stress 
at different points in the structure is expected to occur at 
different times. Some guide as to the probable time of 
occurrence of the most severe stress may be had, however, 
if the computed wing deflection is observed. It is likely 
that maximum stress occurs in the range where wing bending 
appears to be most pronounced. 

The acceleration of any point in the structure may be 
found, if desired, with the aid of equation (5). 

COMPUTATIONAL PROCEDURE 


TABLE 1.— ILLUSTRATION OF THE [S] MATRICES 



The principal results of the analysis presented in the 
previous sections are summarized herein in a step-by-step 
form. Only those steps are listed which actually have to be 
performed when a determination of structural response for 
any airplane is being made. In order to conform with 
standard aircraft practice the use of inch-pound-second 
units throughout is recommended. 

The steps are as follows: 

Preliminary steps: 

(1) The wing semispan is divided into six sections and a 
station is located at the middle of each section (see fig. 3). 
The sections arc proportioned in any convenient manner so 
that certain stations will fall at concentrated mass locations, 
such as engines or fuel tanks. Station 0 is located near 
where the fuselage intersects the wing and station 5 is located 
near the tip. The properties El, GJ, m , rne, and mP are then 
computed at each station. 

(2) From the El, GJ, and X, values determine the [A] and 
[5] matrices by the method given in appendix C. 

(3) Compute the gust-force values at the successive time 
intervals for both the wing and the tail. (See section 
entitled “Derivation of Gust Forces.”) The functions 
used are taken from equations (21) or (22) for the aspect 
ratios which are nearest to those of the wing and tail, respec- 
tively. A time interval that appears satisfactory is one in 
the neighborhood of 1/30 of the estimated natural period 
of the fundamental, bending mode of the wing. 

The recurrence equation: 

(4) With the quantities determined in steps (1) and (2), 
determine the matrix elements given by equations (A3) to 
(A5) at each of the spanwise stations. 

(5) Compute the fuselage and tail coefficients given by 
equations (A8) to (All). (See definitions of M u M 2 , M 3 , 
Mi, M it and M t given by equations (35).) 

(6) With the use of the coefficients determined in steps (4) 
and (5), set up the following matrices: [D\, [S : ], [«.%], [S 3 ] , 
[.R], [e], and \W\. These matrices are defined in appendix A 
and for the most part are found from simple diagonal 
matrices of the coefficients determined in steps (4) and (5). 
The form, for example, of the [S'] matrices is illustrated in 


2 


1 


2 


1 


are not shown are zero. It may be of interest to explain 
briefly the significance of the various terms that appear in 
the matrix. In order to facilitate the explanation the 
matrix has been partitioned into several submatrices. The 
terms in the upper left-hand box are associated with wing 
bending and the airplane vertical moticn; whereas the terms 
in the lower right-hand box are associated with wing torsion 
and airplane pitching. The terms along the two subdiagonals 
serve to couple together the bending and twisting action. 
The terms in the first row and first column are associated 
with fuselage bending. The omission of certain terms in the 
matrix will lead to the matrix which applies to the more 
simple type of aircraft motion. For example, for the case in 
which only wing bending and vertical motion are to be 
considered, computation of only the terms which make up 
the upper left-hand box is sufficient. 

(7) Determine the reciprocal of the [D] matrix and set up 
the following matrix equation: 


S 


w 


= [D]-'\QL 


(7.8) 


where 


V 5 I n 


1 Q\n = 

[-S’,] 

5 

w 

+ [<sy 

5 

w 

+ [S 3 ] 

5 

w 

+ \R\ 





<p 

n — 1 

<p 

n — 2 

<P 

n — ’i 


n 


in which 


\F\ n =[e]\F\,^+\W\ 


In these equations the matrices containing S, w, and y are 
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matrix |F| takes into account the forces which develop due 
to the “wake effect,” and j L g \ is the gust-force matrix which 
is derived in step (3). Equation (78) is seeu to give the dis- 
placements that occur at time n in terms of the displacements 
which occurred at the times n— 1, n— 2, and n— 3. 

The initial response: 

(8) By use of the matrices given in step (6) and the gust 
forces which apply at n= 1, set up the following matrix 
equation : 

I 5 1 


[[T>l + [Sy+8[S 3 ]] 


w 




(79) 


<P li 


The term |F|i does not appear in this equation because it is 
zero. 

(9) Solve equation (79) for the displacements. Any con- 
venient method, such as the Crout method (see reference 6), 
may be used. The displacements found will be the value 
of displacements that apply at t~t or n = 1 . 

The response : 

(10) The response may now be found by successive appli- 
cation of equation (78). The response at n=l has been 
found in step (9); the response at n=2 is next to be deter- 
mined. The values of the displacements in the n—2 term 
of the response equation are all taken to be zero (initial 
condition), and the values in the n — 3 term are taken as the 
negative of those found in step (9). The gust forces to use 
are those which apply at n = 2. The deflections that apply 
at n=2 are then found by matrix algebra. For convenience 
the column matrix |Q| is evaluated first, and then multiplica- 
tion of this column matrix by the reciprocal of the [D] matrix 
gives the deflections at n = 2. With the newly found deflec- 
tions at 71=2 and the deflections at 7i=l and 7i = 0, the 
deflections at ti = 3 can be found, and so forth. This process 
is continued until the wing bending appears to be the most 
pronounced. 

Wing loading: 

(11) With the deflections known, the value of wing loading 
in bending or in torsion can be computed directly from 
equation (77). The stresses at any point can then be com- 
puted from the wing loading by ordinary static means. 
Since the loading may be computed at any value of time, 
the complete stress history of any point on the structure 
may be computed. 



Figure 8.— Weight distribution and equivalent concentrations for example two-engine 

aircraft. 


TABLE 2. — PHYSICAL CHARACTERISTICS AND UNSTEADY 
LIFT FACTORS FOR EXAMPLE AIRPLANE 


b, in - 560 

Co, in - 154 

lb/in..... 165 

p, lb/ft® 0.0765 

T ]mph_. - - - --- 210 

’ jin./sec --- 3700 

v , in./sec 120 

«,see_ 0.01 

As, half-chords - 0. 48052 

M - - 0.276 

At - 10 

tha - - 0.861 

0.001147 

X 0.381 

T- ----------------- - - - 18-3078 

0.832703 

<i»o=ai - 0.361 

<£>„_ -6. 60912 

£ 0 120.9983 


TABLE Z.—i ORDINATES AND GUST-FORCE MATRIX FOR 
EXAMPLE AIRPLANE 


n 

(Equation 

(22)) 

n 

(Equation 

(22)) 

0 

0 

10 

0. 72820 

1 

. 22105 

11 

. 74595 

2 

. 36744 

12 

. 76215 

3 

. 46716 

13 

. 77707 

4 

.53741 

14 

.79086 

5 

.58888 

15 

. 80373 

6 

.62830 

16 

. 81574 

7 

. 65980 

17 

.82696 

8 

. 68594 

18 

. 83748 

9 

.70840 




EXAMPLE 

As an illustration of the method of analysis given in the 
present report, the response of a typical two-engine airplane 
due to a sharp-edge gust is determined. For brevity the 
fuselage is assumed rigid and only vertical displacement 
and wing bending are considered. The weight variation 
over the wing semispan and the equivalent-weight concen- 
trations are shown in figure 8. In this figure are shown also 
the station locations and the interval covered by each station 
section. The solution is made for a forward velocity of 
flig ht, of about 210 miles per hour and a gust velocity of 10 
bet oer second. In tables 2,3, and 4 are listed, respect ively, 


17. 8404 
15. 7552 
12. 1811 

| Lo\ =120 

10. 5295 
8. 77455 
7. 01964 

come from the unsteady lift function, the values of the 
function and the gust-force matrix, and the matrix elements 
that are required for the solution (steps (1) to (5)). The 4> 
function for an aspect ratio of 6 was chosen (see equation 
(20b)); and the function for an aspect ratio of infinity 
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The [A] matrix as computed from the values of X and I 
listed in table 4 is shown in table 5 (a). In the computation 
of the 7] values shown in table 4 for station 0, the fuselage 
was treated as a concentrated wing mass. This treatment 
is allowable since the fuselage is assumed rigid and further 
saves the work of computing the y values (see equations 
(A8)). The [[A] — [iSJ] or [D\ matrix, which in this case 
applies only to bending and vertical displacement, is shown 
in table 5 (b). The equation which is formed from equation 
(78) (step (7)) and which involves the reciprocal of [D] 
and the [Si] and [“/(] matrices is shown in table 6. The 
equation for computing the initial response (step (8)) is 
shown in table 7. 


The solution to these equations is shown in figure 9 in 
which deflection in inches is plotted against spanwise station 
points for various intervals of time. For clarity the deflec- 
tions for the odd intervals have been left off. From these 
curves the consequent wing bending and the manner in 
which the airplane is swept upward by the gust can be seen. 
The time histories of the loads (equation (77)) that occur 
at each of the spanwise stations are shown in figure 10. 
These curves indicate the presence of some second-mode 
excitation in the response. The stresses that occur at stations 
0, 1, and 2 are shown in figure 11. The presence of second- 
mode excitation is not readily discernible from the stress 
curves. 


TABLE 4.— MATRIX ELEMENTS FOR EXAMPLE AIRPLANE 


Station 

X 

I 

c 

l 

Vo 

Vi 

V2 

^3 

S 

0 

0. 09 

2700 

154 

101 

— 56.010712X10* 

139.84200 X10< 

—111. 77100 X10 4 

27.93800 X10 4 

17.9752 

l 

. 18 

1870 

136 

101 

—31. 594032 

78. 802027 

-62. 951014 

15. 733559 

15. 2697 

2 

.17 

1100 

118 

90 

-7. 570015 

18.783512 

-14. 966756 

3. 735946 

12. 2731 

3 

.16 

520 

102 

90 

-2. 109675 

5. 151851 

-4. 060925 

1. 012428 

10. 6091 

4 

. 16 

225 

85 

90 

-1. 150062 

2. 773208 

—2. 168104 

. 539690 1 

8. 84085 

5 

.16 

67.5 

68 

90 

-. 698450 

1.664566 

-1. 291283 

. 320952 

7. 07268 


0 


TABLE 5. — THE [A] AND [D] MATRICES FOR EXAMPLE AIRPLANE 

(a) The [A] matrix 


82, 192. 75 

-133, 410. 07 

61,949. 726 

-12,599.08 

2, 094. 4210 

-237. 6981 

-133, 410. 07 

258, 299. 66 

-172,806. 94 

56, 197. 447 

-9, 339. 8380 

1, 059. 7383 

61, 959. 73 

-172, 806. 94 

194, 219. 495 

—108, 709. 511 

28, 579. 472 

-3, 242. 2365 

-12, 599. 080 

56, 197. 447 

-108, 709. 511 

103, 953. 971 

-47, 410. 326 

8, 567. 4988 

2, 094. 4210 

-9, 339. 8380 

28, 579. 472 

-47, 410. 326 

36, 607. 4681 

-10, 531.197 

—237. 6981 

1, 059. 7383 

-3, 242. 2365 

8, 567. 4988 

-10, 531. 197 

4, 383. 89451 _ 



(b) The [2?) matrix 



642, 389. 82 

-133, 410. 07 

61, 959. 726 

—12, 599. 080 

2, 094. 4210 

-237. 6981 “ 

-133, 410.07 

574, 239. 980 

-172,806.94 

56, 197. 447 

-9, 339. 8380 

1, 059. 7383 

61, 959. 726 

-172, 806. 94 

269, 919. 640 

-108, 709. 511 

28, 579. 472 

-3, 242. 2365 

-12, 599. 080 

56, 197. 447 

—108, 709. 511 

125, 050. 721 

-47, 410. 326 

8, 567. 4988 

2, 094. 4210 

-9, 339. 8380 

28, 579. 472 

-47, 410. 326 

48, 108. 0880 

-10, 531. 197 

-237. 6981 

1,059. 7383 

-3, 242. 2365 

8, 567. 4988 

-10, 531. 197 

11,368.3945 _ 


where 


in which 


TABLE 6.— RECURRENCE EQUATION FOR RESPONSE OF EXAMPLE AIRPLANE 


w(0) 


0. 016448152 

0. 003274436 

-0. 002535913 

-0. 002352327 

-0. 0008201844 

0. 0003284241 ~ 


MJ(1) 


0. 003274436 

0. 022344322 

0. 01467594 

0. 002077166 

-0. 002939533 

—0. 002117350 


w( 2) 

1 

-0. 002535913 

0. 01467594 

0. 069848456 

0. 06246892 

0. 02103792 

-0. 009089948 


ic(3) 

10,000 

-0. 002352327 

0. 002077166 

0. 062468920 

0. 19074974 

0. 15523588 

0. 01762344 


WJ(4) 


-0. 008291844 

-0. 002939533 

0. 02103792 

0. 15523588 

0. 40588417 

0. 26526111 


ra(S) 

n 

_ 0.0003284241 

-0. 002117350 

-0. 009089948 

0. 01762344 

0. 26526111 

1.10968865 _ 






"139. 84200 


w(0) 


"-111. 77100 


tt>(0) 


"27. 938000 


MJ(0) 





78. 802027 


w(l) 


—62. 951014 


w(l) 


15. 733559 


w(X) 


Q 

= 10, 000 


18. 783512 
5. 151851 


w( 2) 
w( 3) 

+ 

—14. 956756 
-4. 060925 


2) 

w(3) 

+ 

3. 735946 
1.012428 


w(2) 

w(3) 





2. 773208 


w( 4) 


—2. 168104 


w(4) 


0. 539690 


w)(4) 



n 


1. 664566_ 


w (5) 

1 

_ -1.291 283 _ 


5) 

n-2 

0. 320952_ 


w(5) 

n-3 


| F\„ =0.832703|F|„ 


"l7. 9752 


w(0) 

15. 2697 


w(l) 

12.2731 


k>(2) 

10. 6091 


W(3) 

8. 84085 


10 (4) 

7. 07268_ 


m>(5) 


TABLE 7.— EQUATION FOR INITIAL RESPONSE OF EXAMPLE AIRPLANE 



" 175. 9720 

-13.34101 

6. 195973 

-1. 259908 

0. 2094421 

-0. 02376981“ 


w (0) 


17. 8404 


-13. 34101 

120. 3415 

-17. 28069 

5. 619745 

-0. 9339838 

0. 1059738 


«*(1) 


15. 7552 

10,000 

6. 195973 

-17. 28069 

41. 92278 

-10. 87095 

2. 857947 

-0. 3242237 


w(2) 

=120(0. 22105) 

12. 1811 

-1. 259908 

5. 619745 

-10.87095 

16. 54357 

-4. 741033 

0. 8567499 


W&) 

10. 5295 


0. 2094421 

-0. 9339838 

2. 857947 

-4. 741033 

6.960225 

-1. 053120 


ui(4) 


8. 77455 


_ —0.02376981 

0. 1059738 

-0. 3242237 

0. 8567499 

-1.053120 

2.413172 _ 


w(6) 

i 

7.01964 


0358 62 -51 3 
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Figure 9. — Response of example airplane due to 10-foot-per-second sharp-edge gust. 
17=210 miles per hour. 




Figure ll. —Bending stress developed in example airplane due to 10-foot-per-second 
sharp-edge gust. 


DISCUSSION 

A method for computing the stresses and structural action 
of an airplane flying through a gust has been given. The 
method is based on aerodynamic strip theory, but over-all 
corrections for compressibility and three-dimensional effects 
can be made as is indicated by a suggested correction pro- 
cedure. Some tail forces are included in the analysis and 
others might equally well be included if their character is 
known. 

The analysis as given is general enough to include the wing 
bending and twisting flexibilities and the fuselage flexibility. 
In a good many cases that may be considered, however, the 
last two of the flexibilities may prove to be. of negligible im- 
portance. Some investigators have indicated (see reference 
1) that, unless the forward speed of the airplane approaches 
the flutter or divergence speed of the wing, the torsional defor- 
mations do not have to be included. Likewise, in cases in 
which the fuselage is rather stiff, the effect of fuselage flexi- 
bility on the response may be rather small. In such cases 
in which either or both of these flexibilities may be ignored, 
the analysis is, of course, simplified and shortened. The 
example presented in the previous section illustrates this 
point. In the present state of understanding of gust- 
response analysis, enough information is not available to 
indicate definitely when and when not to include the various 
flexibilities of the aircraft structure. The analysis in the 
present report may provide a useful means to assess their 
importance. The extent, for example, to which coupling 
exists between wing bending and wing torsion in any par- 
ticular case may be seen by comparing the displacements 
obtained from the coupling terms with the displacements 
obtained from the noncoupling terms. 

Both the symmetrical and antisymmetrical types of gusts 
can be handled by the analysis given in the present report. 
In general, the symmetrical gust is expected to produce the 
most severe stress condition, and therefore only the matrices 
which apply for a symmetrical case have been given. These 
matrices were derived by using the boundary conditions for 
the symmetrical deformation of a free-free beam. The case 
of an antisymmetrical gust can be treated by replacing these 
matrices by the ones which apply for the antisymmetrical 
deformation of a free-free beam. The case of a general un- 
symmetrical gust can be handled by first breaking the gust 
into two parts — a symmetrical part and an antisymmetrical 
part — and then treating each part independently. 

It might be of interest at this point to compare briefly the 
matrix method to a modal-function solution. One of the 
chief disadvantages of the modal-function solution is that the 
modes and frequencies of natural vibration of the structure 
have to be computed in advance. Then, a large number of 
integrals which involve these modes have to be determined in 
order to set up the problem. In the present analysis the 
physical properties of the airplane are used directly in the 
setting-up of the problem. Further, in order to make the 
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modal solution practical the higher modes. must be dropped 
and only the basic or fundamental modes can be used. 
Hence, the success of the analysis depends to a large degree 
on how well single modal functions, one mode each for bend- 
ing- and torsion, can represent the airplane distortion. In 
the analysis of the present report the distortions are found for 
all practical purposes as the correct values at a number of 
spanwise stations, at least to within-the accuracy to whieh-the - 
aerodynamic and structural parameters are known. Also, in 
this analysis, probably the most difficult operation is the 
inversion of the matrix [D\, which is actually not a very 
involved operation, especially when done by the quick and 
systematic procedure afforded by the Crout method (refer- 
ence 6). 

The present report indicates the methods for determining 
the response for both a sharp-edge gust and a gust of arbitrary 
shape. Probably the best approach, however, is to compute 
only the response for a sharp-edge gust, since the response for 
any arbitrary gust may thereafter be computed by means of 
Duhamel’s integral. To follow such a procedure would also 
save a great amount of work in the evaluation of the gust 
forces. 

One of the important features of the method of analysis 
presented is that it is not restricted to the gust problem. 
The approach used may be used to treat other problems of a 
similar nature. The landing problem can be handled by 
simply replacing the distributed gust force by the con- 
centrated landing forces. In the landing problem also, the 
problem is set up much more easily since the aerodynamic 
terms do not ordinarily have to be included. However, the 


landing problem in which aerodynamic forces are included 
may be solved by this method if desired. The approach used 
herein may also be used to solve the problem of the release of 
heavy objects such as bombs. This problem could be consid- 
ered the inverse of the gust problem; a load is released rather 
than-encountered. Some maneuvering problems, such as the 
sudden deflection of the ailerons, and a number of other 
transient problems might also be treated by an approach 
similar to that given in the present report. 

CONCLUDING REMARKS 

A method'for computing, the stresses and structural response 
of an aircraft flying through a gust has been presented. The 
method is based on aerodynamic strip theory, but compres- 
sibility and three-dimensional effects can be taken into 
account approximately by means of over-all corrections. The 
method takes into account wing bending and twisting 
deformations, fuselage deflection, vertical and pitching 
motion of the airplane, and some tail forces. A sharp-edge 
gust or a gust of arbitrary shape in the spanwise or flight 
directions may be treated. A suggested computational pro- 
cedure is given to aid in the application of the method to any 
specific case. 

The possibilities of applying the method to a variety of 
transient aircraft problems, such as landing, are brought out. 


Langley Aeronautical Laboratory, 

National Advisory Committee for Aeronautics, 
Langley Field, Va., January 19, 1950. 


APPENDIX A 

DEFINITIONS OF MATKICES USED IN ANALYSIS 


For convenience in presentation, most of the matrices and 
matrix elements that are used in the analysis are defined in 
this appendix. The matrices are presented without deriva- 
tions, but their origin should become evident by a study of 
the analysis. 

Matrices.— The various matrices that are used in the 
analysis are defined as follows for the case in which the wing 
semispan is divided into six sections (the elements which 
are used in the matrices being defined in the subsequent 
section) : 


M = 


u?(0) 
w(l) 
w(2)\ 
w(3) 
w(4) 
w( 5) 


M = 


*( 1 ) 

<p( 2) 

<p( 3) 
^(4) 
<p(5) 
<p{ 6) 


5 


5 

w 

= 

M 

<p 


bl 


0 


|P| = 


\p\ 

Iffl 


[A] 

[B] 


| See appendix C for definitions. 


[C] = 


“0 0 0 “ 
0 [A] 0 
_0 0 [B]_ 


[V] =[[<?] -[£„]] 


r ’?f(o)+7 i 




Vi( 1) 

Vt(2) 

Vi(3) 

Vi A) 

Vi (5)_ 


~vt(0)+y i 

v/(l) 


lv/] = 


m' (2) 

Vi( 3) 


vi A) 

Vi'(5)_ 


bl = 


p(0) 

pA) 
p( 2) 
P( 3) 
P( 4) 
p{5) 


Vi( 0) + ^< 

F«(l) 


["i] = 


Fi(2) 

v t (3) 


v t A) 


Vi(5) J 


2 ( 0 ) 

2 ( 1 ) 


wm+M/ 

v/(l) 
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■«e 


7 t 

0 



0 


0 


Mi 

0 


! Mi I ~ 


0 

0 


0 






r i o 


L'-iJ^i o o o o oj 
bi'J = Lr/ 0 0 0 0 OJ 
“ ~'i L'’iJL' , i , J _ 
[S ( ]= I "Yi | h i] hi'] 

_ I Mi | hi] hi'] _ 

1 


\m= 


mi m 
i-^i [«-J] 


e~y t t 


e -- ’" 


[e)=\ 


e~y* 


e—>‘ 


0 

0 

0 

0 

0 


— x, 


0 

0 

0 

0 

0 


n 


1 


[/]= 


1 


1 


1 


1 


e _ >' 


[j\=lj ooooo] 

\ j'\=[j’ ooooo] 

-g( 0) 


5 (-l) 


[g}=\ 


g( 2) 


g( 3) 


g(4) 


g(5)J 


k']= 


Y( 0 ) 


s'O) 


g'C 2) 


g'( 3) 


,?'(4) 


g'(5) J 


\W\ = 
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Matrix elements. — The matrix elements which appear in 
the matrices defined in the previous section are expressed 
for convenience in terms of the following common factors: 


P=m A irpU 



$o=«i 

4>o= — 7®i f 

$o=T 2 «i > 


(Al) 


in which the last four are associated with the $ function for 
the wing. (See equation (23) .) With these factors. the ele- 
ments that must he computed at each spanwise station are 
as follows: 


d 0 =— (l-$ 0 ) 0 cl + pcl(* o+^oe) 

di=— (1 — $ 0 ) Pci 
€ 

di= ~Ye (1_$o) ffcl 


<*» = J~ (l-'I'o) Pci 

O € 


ii 0 Z7 e -(io+i^e)c(|-f)] 
= (l-*«)0c*f(§-|) 

d 2 '=l(l-* 0 )PcHQ-<l) 

d.'=-£(i-**p*' i (i-z) 


(A 2) 


2 m , 
170 = 1 — r«o 




5 m , , 
Vi— — 2" + f h 


4m . , 
V 2 = ^2~ + “2 


™+ dl '-^Pcn 

v*=^+d>'+£- e PcH 


y (A 3) 


m . , 

’?3 = ^2 + «3 


v,' = -™+d,'~PcH 



The coefficients which must be computed for the fuselage 
and tail are expressed in part in terms of the following 
common factors: 


Pt = w rn At * P US t 


**o= a ‘« 


2L \ 

yt=~r- 

c "t 


4>, =-7 tdi t 


(A0) 


$i n =7( ®i, 


in which the last four are associated with the <t> functions 
appropriate to the tail. Also used are the generalized 
masses given by equations (35) and the value 81 as given in 
equation (14). With these factors the coefficients for the 
tail and fuselage are as follows: 


/o=-^(i-W+0< (i-o+InO 

/l=f (1 -*«o)0« 

f>=-i v-w* 


(A 7 a) 
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/o'=~ (i-*o (*,+!*)&+* 

^ /= ^( 1 - $, «) (*«+£ C <) ^+ 2 i 7 < J,C ‘ 

y.'=-^ (i-*o (*‘+| c ‘) ^ ,Ci 

9 

/ 2 =-|- 6 ( l -^ o )( 1 +| c < 0 i)^‘- 2 i;i 3 « c ^ 

73=^(i-^ o )( 1 +^ c ' e ‘)^ <+ i^ /3,c<01 


2M, ' 

7o = 72 + /o 


5 M, , r 

Tl=— 2 “+/l 


4M, , j. 

72= T- + /2 


Mi, ' 

73 = - e Y+J 3 

f _ 2M 2 , , 

7o 2 +Jo 


7 i 


5 M 2 . , , 

- 2 - -f-ji 


.. , _ 4 M 2 , j. , 

72 — +J2 


^2 , JT , 

73 — Y+j3 


2 M 3 7 

7o= ~ 2 — I” Jo 

— 5 M 3 | 7 

7i=— 2-+/1 


72 = 


4M 3 7 

,2 rj 2 


— M 3 . 7 

73= — +J3 


(A7b) 


(A7c) 


(A8a) 


(A8b) 


(A8c) 
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Mo = ^V 2 -Z//o 


5 M, 


-xji 


Mi = 


_4M 2 , 

M2 ^2 X tj2 
M3=-^ — Z</3 


Mo 


2M 4 




u' - 5M < x/' 

Ml 2 • c t7l 


M2 


4M 4 




„ ' MA — x f > 

M3 2 


- 2 M’s 

Mo — T 


-a; Jo 


Ml" 


5M 3 , 

2 X Ul 


- 4M 5 - 

M2 — 2 a:(j2 


- _ M 5 7 

M3 2 X V3 


(A9a) 


* 

(A9b) 


(A9c) 


,__2M3 +/o 

^=-^Y 3 +/ 2 

r 3 =^+/ 3 


^o' = -^ 5 +/o' 

5M 5 , , 

? i ~2 I -/l 

/ 4Ms . , , 

' 2 —~~2~-r.T2 


r,' — ¥?+/.' 


2M 6 j !t/ 

2 — h/o— M 6 


7 1 =-f- 6 +/ l 

e" 

7 2 =_iM- 6 + / 2 


~ 6 +/ 3 


j = /3,4>/ 0 te T ‘ e 

1/4>« 0 — ¥, 0 (* ( +| c,)] 
J = Pt ee~ y ‘ e |V4> (o 0i + ¥, o (l+i 0 , 0 ,)] 


(A 10a) 


(A 10b) 


(A 10 c) 


(All) 



APPENDIX B 


DERIVATION OF DIFFERENCE EQUATIONS 


In this appendix the parabolic and cubic difference 
equations for the first and second derivatives of a function 
are derived. 

Parabolic equations. — -For the parabolic difference equa- 
tion, consider the function shown in figure 12(a). This 
function is assumed to be replaced by the arc of a parabola 
which passes through the three ordinates a, b, and c. It 
can be verified readily that such a curve can be given by 
the equation 


d 2 w~\ _ 2d — 5c + 46— a ' 

d y 2 _L = 3c e 2 


(B6) 


If taken at the third of the four ordinates, the derivatives are 


dw~\ _2d + 3c — 6b-\-a 

dy 6« 

d 2 w~ 1 d — 2c +6 

d y 2 _ji/ = 2e e 2 


(B7) 

(B8) 


The first and second derivatives of this equation at y=e are 
given by the equations 


dw~[ _c — a 
dy\,-~ 2 e 

d 2 w ~ 1 c~2b + a 

dy 2 _ L= e — « 2 


(B2) 

(B3) 


These equations are the standard difference equations for 
the first and second derivatives of a function. The deriva- 
tives are purposely taken at the middle of the three ordinates 
because the resulting equations are suitable for use in the 
simplification of many problems. If the derivative had 
been taken at an outer ordinate, the approximation afforded 
would not be accurate enough to be useful. 

Cubic equations. — The cubic difference equations may be 
derived in a manner similar to that for the parabolic equa- 
tions. In this case four successive ordinates are used. (See 
fig. 12(b).) The function is replaced by a third-degree 
curve which is given by the equation 



cy 

2 € 


(!- 0 (!- 2 ) (?- 3 ) + ! t (!~ 2 ) (!- 



(B4) 


Because of the increase in accuracy that results from the 
use of a higher-degree curve, the first and second derivatives 
may be taken at an outer ordinate with an accuracy which 
is about equivalent to that given by equations (B2) and 
(B3). The derivatives at y—3e are 


dw ~ _ 1 Id— 18c + 96 — 2 a 

dyX=Z‘ 


Equations (B5) and (B6) are used in the derivation of the 
response equation for an airplane in a gust. Equations (B7) 
and (B8) are useful in the derivation of the initial response. 




(a) Parabolic. 

(b) Cubic. 

Figure 12. — Functional notation used in the derivation of parabolic and cubic difference 

equations. 
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APPENDIX C 

DERIVATION OF MATRIX EQUATIONS OF EQUILIBRIUM 


In this appendix the matrix equations 

[A] M = |j>| (Cl) 

[B] \<p\ = W | (C2) 


for symmetrical bending and twisting of a free-free beam 
under normal and torsional loads are derived. 

Bending. — In accordance with the assumptions made in 
this report the wing semispan is considered to be divided into 
six sections with a station point at the center of each section 
(see fig. 3). The inertia force of the mass and the aerody- 
namic force that develops over each section is in turn assumed 
to be concentrated at the respective station points. The 
wing is thus effectively a beam bending under six concen- 
trated loads and, as such, will have a linearly varying mo- 
ment between each station. The following general equation 
for the moment between the i and i+ 1 stations may therefore 
be written: 

M=ai + b,y (C3) 

where 




yd) 

b\i + 1 


M{i+ 1) 


The deflection may be found most conveniently from this 
equation by use of the engineering beam theorem which 
states that the deflection of one point on a beam relative to 
the tangent of the deflection curve at another point is equal 
to the moment about the displaced point of the M/EI 
diagram between the two points. In this case symmetrical 
loading is being considered and therefore the boundary con- 
dition at the center line is that the slope must be zero; the 
deflection of each station relative to this point therefore may 
be readily computed. Fortunately, because of the conven- 
ient analytical representation of MjEI, these deflections may 
be found by exact integration. The deflection, for example, 
at station 4 due to the MjEI variation between stations i 
and i+1 may be given by the expression: 

(0i + &.-2/)(c<+d<2/)[i/(4)— 

y(i) 

Consideration of all the expressions of this sort leads to the 
total deflection of each station relative to the wing center 
line. From this deflection the more useful deflection relative 
to station 0 can be readily determined. The values of the 
deflection thus obtained are found to be expressible by the 
following matrix equation: 


o Aj + 1 

in which y(i) is the abscissa to the i station. 

The wing is further assumed to have a linear 1/EI varia- 
tion between stations with the correct value of l/EI at each 
station. This type of variation would lead to an El curve 
which follows very closely the true stiffness curve of the 
wing and which of course has the correct values of El at 
each station. A general equation for \/EI may therefore 
also be written; thus, 


El 


=Ci+d t y 


(C4) 


where 


c*= 1 


L + b\ i+l _\EI(i) 


yii) 


i 


dr- 


EI(i) b\i+iEI(i-\-l) 

l_ r 1 L-i 


b\ i+1 iEI(i+l) EI(i)_ 

With equation (C3) and equation (C4) the well-known 
' expression relating deflection to moment for a beam may be 
written 


^ 2 == Jfj = {<h + b i y)(c i + d { y) 


(C5) 


w(l) 


011 

012 

0 

0 

0 “ 


M( 0) 

w(2) 


021 

0 22 

023 

0 

0 


M{ 1) 

w(3) 

b 2 

031 

032 

033 

®34 

0 


M(2) 

£7(0) 


w(4) 


041 

042 

043 

®44 

045 


M( 3) 

w{ 5) 


_a 5 i 

052 

053 

054 

055_ 


M( 4) 


where the matrix elements are defined by the equations: 

= koki+ Xi 2 /li 

(Z21 = X 0 (Ai + Xo) T- Xi 2 -Ai + X]X 2 I?i 

® 3 i = Xo(Xi-t-X 2 +X 3 )-l-Xi 2 j 4 I -|-X 1 (X 2 -(-X 3 )Bi >- 

1 = X 0 (X i ~t~ X 2 + X 3 -f- X 4 ) + X i 2 Ai + X i (X 2 + X 3 + \j)Bi 
ot 5l = X 0 (X i -i _ A-2 — I - X 3 — f- X 4 — f- X 5 ) -f- X i 2 ^li -)- Xj(X 2 + X 3 — |- X 4 -l— A s )^i ^ 

(C7a) 

012=Xl 2 Ci 

02 3 == X i^C\ -f- X iX 2 Hi -j- X 2 2 A-2 

032 = X i 2 Cj — f-XifXj— f-X 3 )Z?i-f- X 2 2 A 2 TX2X 3 I?2 

042= X i 2 Ci — F X i(X 2 — |— X 3 — X 4 ) Z?i — |— X2 2 ^A 2 -f— X 2 (X 3 — f— X 4 ) 

052^= Xi 2 <Ti-(— \i(x 2 -(— x 3 — (— x 4 -f- x 5 ) x 2 2 -^l 2 + x 2 (x 3 — (-x 4 -f- x 5 ) a? 2 _ 

(C7b) 
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in which 


G&23 ^2 2 f') 

a 33 =\ 2 *C 2 + \ 2 -\ 3 D 2 +\ 3 *A 3 

Cr-43 = X2 2 C I 2"1 - ^2(^3 "t - X 4 ) D% X3 2 A3 -f- X3X4-B3 

a 53 ~ X 2 2 C 3 -f" X 2 (X 3 -j- X 4 + X 5 ) Z ?2 "I" X 3 2 44.3 + Xs(X 4 + X 5 )i ?3 ^ 

(I44 = \^C 3 -f" X3X4ZI3 X4 2 -44 

<154= X3 2 C 3 + X 3 (X 4 + X5) Z) 3 + X 4 2 yl4 + X4X5 B 4 ^ 

GE'45=X 4 2 0 4 "1 

, ^55 ~~ X 4 2 C \ ~ (“ X 4 X 5 2^)4 " 4 " A . 5 J 

1 1(0) 1 7(0) 'I 

' 4 I(i — 1) ‘ 12 I(i) 

R 1 1(0) | 1 HO) 

1 3 /(i— l)~o I(i) 

r 1 1(0) 1 /( 0 ) ” 

'* 12 I(i — 1) ' 12 I(i) 

1 7(0) , 1 7(0) 

1 6 I(i — 1) ‘ 3 I(i) 


(C7c) 


(C7d) 


(C7e) 


(C8) 


The moment M(5) does not appear in equation (C6) because the boundary condition is used so that the moment at station 
5 is zero. For convenience equation (C6) may be given in the abbreviated form: 


\™\=Em) [//i] |M| 


(C9) 


From the five loads p( 1), p( 2), p( 3), p( 4), and p( 5), the moment at each station may be found. The equations 
relating the moments to the loads can be shown to be given by the matrix equation: 


M(0) 


"X! 

X1 + X2 

X1 + X2+X3 

Xi + X 2 +X 3 +X 4 

Xi + X 2 + X s + X 4 + X 5 


p(l) 

M{ 1) 


0 

x 2 

X 2 + X3 

x 2 + X3T X 4 

X2 + X3 + X4+X5 


p( 2) 

M{ 2) 

= b 

0 

0 


X3+X4 

X3+X4+X5 


m 

M(3) 


0 

0 

0 

X4 

X4+X5 


p( 4) 

M( 4) 


_0 

0 

0 

0 

Xs_ 


PC 5) 


(CIO) 


which can be abbreviated simply 


results in the equation: 


\M\=bm\ P \ 


(Cll) 


Substitution of equation (Cll) into equation (C9) gives 




(Cl 3) 


\”\=EI( o) ^ ™ \P\ 


(Cl 2) 


Multiplication through by the reciprocal of 


EI{ 0) 


[//.] [HA 


This equation thus gives the loads in terms of the deflection 
of each station relative to station 0. In the case under 
consideration, however, the wing is a free body capable of 
motion through space and therefore to set up properly the 
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equations of motion the deflection must be given relative 
to a fixed datum line. This datum line is most conveniently 
located as the position of the wing prior to action of the 
disturbing loads. Consideration of sketch 3 


£ 



will show therefore that the following relation must exist: 

w=w—w{ 0) (Cl4) 


Equation (C17) is noted to express all the loads except 
p(0) in terms of the six deflections. An additional equation 
in which p( 0) is expressed also in terms of the deflections 
may be established by use of the condition that all the loads 
acting on the wing semispan must add up to equal zero; 
that is, 

7(0)+z( 1 ) + Z(2)+z(3)+z(4) + z(5) = 0 (Cl 9) 

This condition automatically satisfies the two boundary 
conditions that the shear must be zero at the tip and center 
line of the wing. Thus if the five equations represented by 
equation (Cl7) are added, and use is made of equations 
(Cl 8) and (Cl 9), the following equation results: 

&oow(0) + 6 0 ,w(l) + 6 02 w(2) + 6 0 sw(3)+ 6 04 w(4) + b 05 w( 5 ) = p(0) 

(C20) 

where 

boo— — (6qi + 6 02 + 603 + 604 + 605) (C21) 


Substitution of this equation into equation (Cl3) gives 

[[#,] [# 2 ]]->— w( 0 )| = bl (C 15 ) 

j? T((\\ 

To aid in the derivation — —■ [[//,] [Z^]] -1 is now written in 
the general form: 

bn b I2 b i3 b i4 6 15 

b 12 b 22 b 23 b 24 b 25 

[ 771 [Z 7 ]] I= | b 13 623 633 634 635 j (Cl 6) 

614 b 24 b 34 b 44 b 45 

b 15 6 2 5 b 35 b 45 655 

Thus with this equation, equation (Cl 5) may be trans- 
formed to the form: 


£7(0) 


where 








w { 0) 


boi 

bn 

b 12 

b 13 

6 14 

6,5 


w(l) 


7 ( 1 ) 

b 02 

b 12 

b 22 

b 23 

6 2 4 

625 


w { 2) 


Pi 2 ) 

b 03 

b 13 

b 23 

b 33 

bu 

635 


w ( 3 ) 

= 

pi 3 ) 

b 04 

bu 

624 

bn 

644 

645 


w ( 4 ) 


pi 4 ) 

b 05 

bin 

6 2 5 

b Ob 

bib 

655 


w { 5 ) 


Pi. 5 ) 


b oi= 

~(b 

11+6 

12+6 

13+6 

14+615) 



(Cl 7) 


b 02= — (612+ 6 2 2 + 6 2 3 + ^ 24 1 625) 

b 03= ((>13 "I - b 23-)- 633+ 634+ 6 35) 

f>04— '~(6l4 + 6 24 + ^ 34 I 644 + (*45) 
boi — — ((>15+ 6 25 + 635+645 + 655) 


(Cl 8) 


This equation may now be combined with equation (Cl7) 
to give finally 


boo 

601 

602 

603 

604 

605 


w(0) 


7(0) 

601 

6„ 

612 

613 

6 14 

615 


w{l) 


7(1) 

602 

612 

6 22 

623 

624 

625 


w{2) 


7(2) 

603 

613 

.623 

633 

634 

635 


w ( 3) 


Pi 3) 

604 

6,4 

624 

634 

644 

645 


w(4) 


7(4) 

_ 605 

6,5 

625 

635 

645 

6 5 5_ 


w(5) 


7(5) 


(C22) 


This equation is thus the desired matrix equation which 
relates the normal loads to the deflection. If the square 
matrix is denoted by [A], the equation may be abbreviated 
conveniently to the form 


[A ] M = |zl 


(C23) 


which is the form used in the text. (See equation (40).) 

As an aid in computational work, a summary of the steps 
involved in the determination of [A] is given to close this 
section: 

(1) From the I values at the respective stations, compute 
the coefficients given by equations (C8). 

(2) With these coefficients determine the matrix elements 
given by equations (C7). These elements form the matrix 
[17] which is defined by equations (C6) and (C9). 

(3) Multiply the [77,] matrix by the 77] matrix, which is 
defined by equations (CIO) and Cll). The result should be 
a symmetrical matrix; this property serves as a very useful 
computational check. 

b 3 

(4) Invert the [ 77] 77] ] matrix. This matrix should 

also be symmetrical. (The Crout method (reference 6) 
serves as a rather quick and useful means for performing the 
inversion.) 
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(5) Add the columns of the inverted matrix and place the 
negative of these sums at the top of their respective columns 
such as to form a new row of matrix elements. Then add 
these sums and place the negative 5f the sum as the first 
matrix element of the newly formed row. A new column 
headed by this value is thus in the making. Fill in the re- 
mainder of the column with the respective elements of the 
new row; that is, the appropriate values should be inserted to 
make the matrix symmetrical. This final matrix is the de- 
sired [A] matrix. 

Torsion. — For the torsional case the torque loads q are 
assumed to be concentrated at the stations just as in the 
case for the normal loads p. Consideration then of the fol- 
lowing example torque diagram (sketch 4) 

<L 



1 

1 

\T(Z) 

1 T ( 3 / 

T(4) 

1 T(5) 

, ^ ^ r. / 


| V. — c V- X3 


9(0) q(V 9(0) q(3) q(4) q(3) 

Sketch 4. 

will show that the following equations must apply: 

2 (0)= -T{ 1) 

q{l) = ni)-T{2) 

2(2) — T(2) — T(S) 

2(3)=T(3)-T(4) 

2(4)=T(4)-T(5) 

2(5)= T(5) 

where T(i) represents the total torque present in the i interval. 
No torque exists between the wing center line and station 0. 

To aid in the derivation, the assumption is made that 
1/GJ varies linearly between stations. A typical TjGJ 
diagram between, say, the i — 1 and the i stations would 
appear as in sketch 5: 


y 

(C24) 


d T 

From the differential relation the feet may be ob- 

served that the change in angle of twist between two stations 
is equal to the area of the T/GJ diagram between the two 
stations; therefore, 




2 lGJ(i-l) 


If the notation 


. 2 G 1 

3 ‘ b\t 1 1 

J(i-1) + J(i) 


m i 

GJ{i)\ 


(C25) 


(C26) 


is employed, equation (C25) may be written 


T{i)=j t {<p{i)-<p{i- 1)] (C27) 


Application of this equation to each of the spanwise stations 
gives the following equations for T: 


T(l)=ji Wl)-^(O)] 
T(2)=j 2 [*(2)-*(l)l 

T(4)=j 4 [«,(4)-«>(3)] 

7 1 (5)=,7 5 [<t>(5) <p(4)] „ 


(C 28 ) 


Substitution now of these equations into equations (C24) 
gives the desired equations relating the torque loads to the 
angle of twist. The equations thus found can be given in 
the matrix form: 


A 

-ji 

0 

0 

0 

0 " 


v(P) 


2(0) 

~ji 

(il+ 7 ’ 2 ) 


0 

0 

0 


<P( 1) 


2(1) 

0 

— )2 

(ji+jj) 

-is 

0 

0 


¥>(2) 


2(2) 

0 

0 

—is 

(is+id 

—ji 

0 


v(3) 


2(3) 

0 

0 

0 

—ji 

(ji+ji) 

— js 


<p(4) 


2(4) 

0 

0 

0 

0 

-is 

j 5_ 


<p(5) 


2© 


(C29) 

which can be abbreviated to 


m 


Tfi) 

GJfi) 

GJ(i-l) 

bXi 



5 5* 



i-l i 


[£]M = l2l (C30; 

the form used in the test. (See equation (41).) Thus all 
that is involved in the computation of the matrix [ B ] is the 
evaluation of the matrix elements by means of equa- 
tion (C26). 


Sketch 5. 



APPENDIX D 


RECURRENCE EQUATION FOR THE EVALUATION OF DUHAMEL’S INTEGRAL INVOLVING AN EXPONENTIAL KERNEL 

9 


The derivation of a rather simple recurrence relation for 
the step-by-step evaluation of the three unsteady lift inte- 
grals appearing in equation (25) is presented. This deriva- 
tion is made possible because the kernels of the integrals are 
expressible in exponential form. 

From equation (23) the first and second derivatives of the 
$ function may be written 


- 217 -x^, ■ 

<f>= \a x e c ° 

6 0 

(Dl) 

.. 4[/2 -X .. 

c » =$ 0 e- yl 

(D2) 


where 


$ 0 = — ya 1 

$o=7 2 «i 

With these equations the three integrals of equation (25) 
may be combined conveniently into the following single 
integral denoted by I,: 


where F 0 does not appear since the initial conditions are 
used that the deflection w and rotation <p are zero at f— 0, 
and therefore Fo is zero. (See equation (D4).) More ac- 
curate methods, such as Simpson’s method, could be used 
for determining the area under the curve, but because of the 
small interval chosen the consequent increase in accuracy is 
negligible. If the notation 

f, = t 6- w [7 1 e" + 7 ! e*-l + (D7) 

is introduced, equation (D6) may be written simply 

In=F n +±Y n (D8) 

If equation (D5) is expanded similarly, only for an upper 
limit of t— e, the expanded result would be 

=ee- y(n ~ lu F 2 e T2< 4 + Y n ~ 2 e y(n - 2)e + 

\ (D9) 

By analogy with equation (D7), however, 
F„_ 1 =ee- T(B - 1)6 [F 1 e*+F 2 e 1 ' 2e - f • • • + F B _ 2 e T <"- 2 > e ] (DlO) 


I l= J o faclw-~[*oPclU+* 0 (3cH Q-j)] „ jc * (l 

(D3) 

For convenience the notation 

F = ^oPclw—^QnPdU +$ 0 Pc 2 l (| — f )J ^ j (D4) 

is introduced and thus equation (D3) becomes 

I,=j'Ye-y (t -*d T 

or 

= Ye y, dT (D5) 


and therefore equation (D9) becomes 

I»-\= Fn-l-\~2 r„_i (Dll) 

A study of equations (D7) and (DlO) shows that the follow- 
ing relation must exist: 

F n =e-r*F n _ 1 +ee-y e Y a - l (Dl2) 

Now, if equation (D4) is used to rewrite Y n and Y n _ x in 
equations (D8) and (Dl2), the value of /„ may be given 
finally by the equation: 

In = F n + i $ 0 |8c( ew„ — ~ Pel e C 0— f <Pn (D 1 3 ) 


Mathematically, the integral in this equation may be 
interpreted to represent the area under the function given 
as a product of F and e yT . .In accordance with numerical 
evaluation processes, the interval 0 to t may be divided into 
a number of time stations of interval e. The product of 
F and e yT may then be found at each of the time stations 
and from these products the area under the curve may be 
determined in first approximation by the trapezoidal method 
of determining areas. Thus, if the time station n corresponds 
to time t, the expression for /, may be approximated as 
follows: 

I t ^I n =te- yn '^Y l e yt -\-Y 2 e y *+ . . . + Y n _ 1 e y< - n - lu +^Y n e yne J 

(D6) 


where 

F n = e ~ y< F n _ i + $ 0 te ~ y ‘Pclw n _ , — 

(D14) 

The value of the unsteady lift integrals is thus given by 
equation (D13). As regards the analysis given in the 
present report, w n -i and are the values of deflection and 
rotation which have, say, just been determined from the 
recurrence equation for response. The value F n _, was also 
established and therefore F n can be determined as a definite 
quantity. The value I n is thus seen to be given in terms of 
the known F„ and in terms of w n and <p n which are the next 
values to be evaluated from the recurrence equation. 



APPENDIX E 

MATRIX ALGEBRA 


This appendix is written for those not familiar with matrix notations or matrix methods. All the matrix algebra necessary 
for the understanding of this report is described hereinafter by way of examples. 

Matrix definition. — Some of the basic types of matrices are illustrated by the following arbitrary matrices which are of 
the third order: 


The column matrix 


The row matrix 
The square matrix 


The diagonal matrix 


The identity matrix 


2 

1 



-1 


L 2 

-3 

U 

" 2 

—3 

r 

1 

2 

— 2 

1 

— 1 

3 _ 

" 4 

0 

0" 

0 

3 

0 

0 

0 

— 1 _ 

" 1 

0 

0” 

0 

1 

0 

0 

0 

1 _ 


Element definition.— Each of the terms that appear in a matrix is defined as an element. Its position is usually 
denoted in a row by the number of terms from the left and in a column by the number of terms from the top. 

Matrix addition. — The addition of two matrices produces a single matrix. Addition is performed by simply adding 
together corresponding elements. For example, 


" 2 

-3 

1 “ 


'4 

— 1 

0 " 


'6 

— 4 

1 " 

1 

2 

— 2 

+ 

0 

3 

2 

= 

1 

5 

0 

_-l 

— 1 

3 - 


_2 

0 

-1- 


_1 

-1 

2 - 


Multiplication of a matrix by a scalar number. — In the multiplication of a matrix by 
the matrix is multiplied by the number. For example, 


a scalar number every element in 


” 2 

-3 

1 " 


' 4 

-6 

2 " 

1 

2 

— 2 

= 

2 

4 

-4 

1 

— 1 

3 - 


_ — 2 

-2 

6 - 
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Multiplication of a column matrix by a row matrix. — The product of a column matrix and a row matrix is equal to the 
sum of the products of the corresponding elements. For example, 


12 -3 1J 


2 

1 

— 4 j 


= ( 2 X 2 ) + (- 3 Xl) + [lX(- 4 )] = -3 


Multiplication of a column matrix by a square matrix. — The multiplication of a column matrix by a square matrix 
produces a column matrix. Consider the following set of three simultaneous equations: 


2 yi— 3 1/2+ ?/3 =«r 
yi + 2 y 2 — 2 r/3=a 2 f 
—Vi— y2 + 3y 3 =a s ^ 


(El) 


The procedure adopted in matrix algebra is to write these equations in the matrix form 


i-H 

CO 

1 

CM 

1 


2/i 


a, 

1 2 —2 


2/2 

= 

do 

1 —1 3_ 


2/3 


d Z 


(E2) 


where the multiplication of the \y\ column matrix by each row in the square matrix produces the respective elements 
in the |a| column matrix. (See multiplication of a column matrix by a row matrix.) 

In order to simplify the presentation of an analysis, the symbolic or abbreviated matrix form is used quite often. 
The symbolic form of equation (E2) is simply 

[M]\y\ = \a\ (E3) 


The determination of \a\ by the multiplication of \y\ by [M] is illustrated with arbitrary values of y, say 2/1 = 4, 2/2 = 5, 
and 2/3=6, by the following equation: 


r 2 -3 

1 ■ 


4 


(2X4) + (-3X5)+ (1X6) 


— 1 


a, 

1 

2 

—2 


5 

= 

(1X4)+ (2X 5) + ( — 2X6) 

= 

2 

= 

a 2 

_ — 1 

— 1 

3 _ 


6 


(-lX4) + (-lX5)+ (3X6) 


9 


a 3 


Multiplication of a square matrix by a square matrix. — The multiplication of two square matrices produces a square 
matrix. Multiplication is performed by letting the multiplying matrix operate, as in the preceding section, on each 
of the successive columns in the matrix being multiplied to produce corresponding successive columns in the product 
matrix. For example, 


"2—3 1 " 


" 1 —2 3" 


“—2 —12 14" 

1 2 —2 


1 3 —2 

= 

5 2-5 

_-l -1 3 _ 


1 1 2 _ 


_ — 5 2 5 _ 


(E5) 


Order of multiplication. — In general 
matrix methods; that is, 


the commutative multiplication law of ordinary algebra does not hold in 


\a\\b\*\b\ \a\ 


Therefore, whenever the product of several matrices is indicated, these matrices must be multiplied together without 
changing their order. 
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Matrix partitioning and submatrices. — A matrix may be 
partitioned or divided at will into smaller matrices For 
example, the left-hand side of equation (E4) may be parti- 
tioned as follows : 


" 2 j -3 1* 


4 

1 1 2 —2 


5 

1 ! -1 3_ 


6 


The matrices which are formed by the dividing lines are 
called submatrices. These submatrices may be treated 
as though they were elements when matrix operations are 
performed. For example, with the notation 

o=[— 3 ij 


where [M\~ l is the reciprocal, or the inverse, of [M], 

The reciprocal of a matrix is found as the matrix which 
satisfies either of the equivalent equations 

[M]- 1 [M] =[/] 

[M\ [M] _1 =[7] 


where [/] is the identity matrix. For equations (E2) and 
(E3), the reciprocal of [M\ is found as the matrix which 
satisfies the equation 


1 

i-H 

CO 

1 

<N 

1 


b\ C\ d\ 


0 

0 

i- i 

1 

1 2 —2 


62 C2 <^2 

= 

0 1 0 

l — 

1 

1 

H-* 

CO 


_ 63 £3 d%_ 


0 

0 

1 


1 

b = 

-1 



5 

d = 

6 


the multiplication of the foregoing partitioned matrix is as 
follows : 


" 2 j a~ 14 


8 ad 

b j c_ d 


4 b + cd 


If this equation is considered in relation to equations (El), 
(E2), and (E5), the values of b u b 2 , and b 3 would simply be 
values of y u y 2 , and y 3 which satisfy equation (El) for 
a,= l, a 2 = 0, and a 3 =0; C\, c 2 , and c 3 would be the values for 
a 2 = 0, a 2 =l, and a 3 = 0; and d u d 2 , and d 3 would be the 
values for a, = 0, a 2 = 0, and a 3 — 1. For this example, the 
solutions are 


t 1 _ 7 , _ 5 

" 2 12 C2 12 “ 2 12 

i 1 _ 5 ,7 

63 12 Ci 12 ® 3 12 


The reciprocal of a matrix and the identity matrix. — By 

ordinary algebraic methods the formal operation involved 
in the solution for x of the equation 

mx=a 


The Crout method (reference 6) provides a very quick and 
convenient method for determining these solutions. 

The determination of y by the operation [M\~ l on \a\ is 
illustrated as follows for a x — — 1, a 2 = 2, and a 3 = 9: 


, is the multiplication through by the reciprocal of m; thus, 

x—m~ 1 a 

The same formal operation may be applied to matrix equations. 
For example, the solution for \y\ in equation (E3) is simply 

|y| = [MT I |a| 



(Ell) 


The operation performed by this equation can be seen to be 
the inverse operation of equation (E4). 
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Positive directions of axes and angles (forces and moments) are shown by arrows 


Axis 

Force 
(parallel 
to axis) 
symbol 

Moment about axis 

Angle 

Velocities 

Designation 

Sym- 

bol 

Designation 

Sym- 

bol 

Positive 

direction 

Designa- 

tion 

Sym- 

bol 

Linear 
(compo- 
nent along 
axis) 

Angular 

Longitudinal 

X 

X 

Rolling 

L 

Y *Z 

Roll 


u 

P 

Lateral 

Y 

Y 


M 

z — >x 

Pftfth 

0 



Normal 

z 

z 


N 

X >Y 


i b 

H 

f 










Absolute coefficients of moment 

ri _ j* r< Af « N 

v, ~$bS *~%eS H ~gbS 

(rolling) (pitchmg) (yawmg) 


Angle of set of control surface (relative to neutral 
position), 5. (Indicate sin-face by proper subscript.) 


4. PROPELLER SYMBOLS 


D 

Diameter 

T) 

V 

Geometric pitch 

Jr 

P/D 

Pitch ratio 

c. 

V' 

Inflow velocity 

V. 

Slipstream velocity 

V 

T - 

Thrust, absolute coefficient C T = — ttfh 

pnU* 

n 

Q 

Torque, absolute coefficient O, Q~~ p tfjy& 

$ 


Power, absolute coefficient ^p— pn sj)i 

Speed-power coefficient = -y/p~“* 
Efficiency 

Revolutions per second, rps 
Effective helix angle^tan -1 ^^^) 


1 hp=76.04 kg-m/s=550 ft-lb/sec 
1 metric horsepower =0.9863 hp 
1 mph=0.4470 mps 
1 mps=2.2369 mph 


5. NUMERICAL RELATIONS 

1 lb=0.4536 kg 
1 kg=2.2046 lb 
1 mi= 1,609.36 m=5,280 ft 
1 m=3.2808 ft 


L 





